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1.1 
Introduction 



o 

^ Attosecond science has emerged with the discovery of coherent electron-ion coUi- 

■ ^ sions induced by a strong laser field, usually referred to as "re-collisions" (ICorkuml 

C/3 



(|1993)). This discovery was initiated by the numerical experiments of K. Schafer, J. 



(-H Krause and K. Kulander (see|Krause et fl/.|(|1992a[l). The work by Corkum (|1993[l 



O , drew on the concepts developed in the earlier work of |BruneI| ([1987 1990 1 and 



[Corkum et al.\ ([T989]f. It has also been predated by the concept of the 'Atomic 



Antenna' ([Kuchiev ( 1987)). With the benefit of hindsight, we now see the work by 



^ [Kuchiev'(1987"l as the earliest quantum counterpart of the classical picture developed 

CT^ by [Corkum (^993^1 and Kulander et a/.| ( |l993l lp| 



^~~^ The classical picture of strong-field-induced ionization dynamics is summarized 

^^ as follows. Once ionization removes an electron from an atom or a molecule, this 

electron finds itself in the strong oscillating laser field. Newton's equations of mo- 



\l tion show that, within one or few cycles after ionization, the oscillating electron 

^~^ can be driven back by the laser field to re-encounter the parent ion. During this 

I re-encounter, referred to as re-collision, the electron can do many things: scatter 

• • elastically (diflract), scatter inelastically (excitation or ionization of the parent ion), 

. 5^ or radiatively recombine into one of the ion's empty states. It is this latter process 

^^ that we will focus on here. The classical picture is usually referred to as the three-step 

5-H model, or the simple man model rj 

If the recombination occurs to the exact same state that the electron has left from, 

1) While the quantum vision of [Kuchiev ( 1987 1 has predated the classical picture, at that time it lacked 
the striking clarity and transparency of the classical model ^Corkum [ ^ 1 993| ), which linked several key 
- and seemingly disparate strong-field phenomena - high harmonic generation, production of very high 
energy electrons, and extreme efficiency of double ionization. The history of this discovery is rich and 
interesting in its own right. Some of it is recounted, from a more historical perspective, in Chapter 4. 
Our purpose here is different - we simply urge our reader to read the papers by |Brune'l[^1987|[l99Q^ , 
jCorkum et al. 1 1 989) , [K'uchiev|(T987) , [Schafer et aL\[\99'i) , as well as a seemingly uiu'elated paper of 
['Gallaghernl988f . 

2) As far as one of us (M.I.) can remember, the latter term has been used by K. Kulander, K. Schafer and 
H.-G. MuUer, who have contributed a lot to the development of this classical model. 
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then the phase of the emitted radiation is the same from one atom to another, leading 
to the generation of coherent radiation in the medium. This process is known as high 
harmonic generation (HHG). It produces tens of eV-broad coherent spectra and has 
two crucial appUcations. First, high harmonic emission is used to generate attosec- 
ond pulses of light (see e.g. [Krausz and Ivanov|p009[ l), which can then be used in 
time-resolved pump-probe experiments. Second, the ultra-broad coherent harmonic 
spectrum carries attosecond information about the underlying nonlinear response, 
which can be extracted. The second direction is the subject of high harmonic spec- 
troscopy (see e.g. 'Lein'(2005"i;"Bak er era/.| ( (2006) ; |Smirnovaef a/.| ( |2009a^ ; |Haessler| 
\et al. | ( [20T0 j ) - a new imaging technique with a combination of sub- Angstrom spatial 
and attosecond temporal resolution. 

In the language of nonlinear optics, high harmonic generation is a frequency up- 
conversion process that results from the macroscopic response of the medium. The 
nonlinear polarization is induced in the medium by (i) the response of the atoms and 
the molecules, (ii) the response of the free electrons, (iii) the response of the guiding 
medium. Here we focus on the theory of single atom or single molecule response. 
The description of macroscopic propagation effects, which determine how coherent 
radiation from different atoms or molecules add together, can be found in |Gaarde| 
|eFar] ( |2008j . 

From the famous simple man model to the recent multichannel model, we will try 
to guide you through the several landmarks in our understanding of high harmonic 
generation. We hope to provide recipes and insight for modelling the harmonic re- 
sponse in complex systems. The chapter includes the following sections: 

1.2 The simple man model of high harmonic generation (HHG); 

1.3 Formal approach for one-electron systems; 

1.4 The Lewenstein model: stationary phase equations for HHG; 

1.5 Analysis of complex trajectories; 

1.6 Factorization of the HHG dipole: simple man on a complex plane; 

1.7 The photoelectron model of HHG: the improved 'simple man'; 

1.8 The multichannel model of HHG: Tackling multi-electron systems; 

1.9 Outlook; 

1.10 Acknowledgements; 

1.11 Appendix A: Supplementary derivations; 

1.12 Appendix B: The saddle point method; 

1.13 Appendix C: Treating the cut-oflF region: regularization of the divergent sta- 
tionary phase solutions; 

• 1.14 Appendix D: Finding saddle points for the Lewenstein model. 

Atomic units h = m — e — 1 we. used everywhere, unless specified otherwise. 
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1.2 

The simple man model of high harmonic generation (HHG) 

Experiments in the eighties and the early nineties of the last century yielded an as- 
tounding result: shaken with sufficiently intense infrared laser radiation, the atomic 
medium was found to up-convert the frequency of the driving infrared laser light by 
up to two orders of magnitude (see e.g. |Huillier and Balcou| ( |1993[ l; |Macklin et al.\ 
( |1993) ). The observed harmonic spectrum formed a long plateau, with many har- 
monic orders, followed by a sharp cut-off. This observation has to be placed in the 
context of what has been routinely seen in the traditional nonlinear optics: in the 
absence of resonances, the nonlinear response would decrease dramatically with in- 
creasing harmonic order, and the harmonic numbers would hardly ever reach double 
digits, let alone form a plateau extending beyond N=101. 

To generate very high harmonics of the driving frequency, the atom has to ab- 
sorb lots of photons. Generation of harmonics with numbers like N=21,...,31,..., etc. 
means that at least that many photons (21, ..., 31, ...) had to be absorbed by the atom. 

The minimal amount of photons required for ionization is A^q = Ip/u, where Ip is 
the ionization potential and oj is the infrared laser frequency. For Ip ~ 12 — 15 eV and 
an 800 nm driving IR laser field (the standard workhorse in many HHG experiments), 
N{) ~ 10. One would have thought that once ten or so photons are absorbed, the 
electron should be free. And since it is well-known that a free electron should not 
absorb any more photons, the emission should stop around A'^ — 11 or so, in stark 
contrast with experimental observations. 

Why and how many additional photons are absorbed? What is the underlying 
mechanism? The liberated electron oscillates in the laser field, and its instantaneous 
energy can be very high. Can this instantaneous electron energy be converted into the 
harmonic photons? Where is the source of non-linearity, if the free electron oscillates 
with the frequency of the laser field? 

The physical picture that clearly answered these questions is the classical three-step 
model. It is simple, remarkably accurate, and is also intrinsically sub-cycle: within 
one optical period, an electron is (i) removed from an atom or molecule, (ii) accel- 
erated by the oscillating laser field, and (iii) driven back to re-collide with the parent 
ion. This picture connects the key strong-field phenomena: above-threshold ioniza- 
tion, non-sequential double ionization, and high harmonic generation. It reveals the 
source of non-linearity in HHG: the recombination of the accelerated electron with 
the ion. 

How can one check that this mechanism is indeed responsible for HHG? The key 
thing test is whether or not this picture explains the cut-off of the harmonic spectra, 
that is, the highest harmonic order that can be efficiently produced. Numerically, the 
empirical cut-olTlaw was found to be fimax ~ Ip+3Up ( [Krause e? a/. | ( | 1 992b J ), where 
Up is the cycle-averaged energy of the electron's oscillatory motion in the laser field. 
To calculate the classical cut-oflF, we should calculate the maximal instantaneous en- 
ergy of the returning electron, but to do so we need to know the initial conditions for 
the electron just after ionization. These conditions are specified within the three-step 
(simple man) model of HHG, which makes the following assumptions: 
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Figure 1.1 Window of classical 'birth' times and the return energy. Left panel: Time of 
birth vs. time of return. Right panel: Energy of the electron at the time of return. 



• SMI: The electron is born in the continuum at any time within the laser cycle; 

• SM2: The electron is born near the ionic core (i.e., near the origin of the reference 
frame) with zero velocity; 

• SM3: If the electron returns to the ionic core (origin), its instantaneous energy at 
the moment of return is converted into the harmonic photon. 

The pull of the ionic core on the liberated electron is neglected in the model, which 
is not unreasonable considering the very large excursions that the electron makes in 
the strong driving laser field. The possibility of the electrons return to the core is 
dictated by the phase of the laser field at which it is launched on its classical orbit, 
and the time-window for the returning trajectories - the range of the 'birth' times tg 
- shown in Fig. |l.l[ 

The calculation is done as follows: for each t^, we find the time of return tj^ to 
the electron's original position (Fig. |l.l| left panel) and the energy at the moment of 
return (Fig. |l.l| right panel). The assumption that the strong laser field dominates 
the electron's motion after ionization simplifies our calculations. Once the ionic core 
potential is neglected, the kinetic momentum (velocity) at the time of birth tg can be 
written as k(ts) = P + A(fB), where p is the canonical momentum of the electron 
and A(t) is the vector potential of the laser field, which is related to the electric field 
F{t) as F(t) = ~dA/dt. The condition k(tB) = (SM2) specifies p = -A(tB). 
Therefore, the electron kinetic momentum at all later times t is k(i) — —A{tB) + 
A(i) and the electron energy at the time of return is 

Erctitn) = k2(i^)/2 = (A{tB) - A{tn)f/2 . 

The zero displacement of the electron from the time of birth, ts, to the time of return, 
tfi, (SMS) defines the return time f^: 

dt{A{t)~A{tB)) = Q. (1.1) 



According to this model, the maximal return energy is about 3.17 Up, where Up — 
F"^ /AuP' and F is the electric field amplitude (see Fig. 1.1). Then the maximum 
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energy of the emitted harmonic photon is 3.17 C/p+/p, where Ip is the binding energy 
of the ground state to which the electron recombines, is in excellent agreement with 
the empirical cut-off law found numerically by |Krause et a/.|p992a) . 

The formal quantum approach considered in the next section will first take us away 
from the simple classical model. However, just like the re-colliding electron revisits 
the ion, we will revisit the simple man model several times in this chapter, refining it 
at each step. 



1.3 

Formal approach for one-electron systems 

The response of an individual atom or a molecule P (r, t) — nD (i) is proportional to 
the induced dipole D(t): 

D(t) = (*(t)|dl*(i)), (1.2) 

where n is the number density, d is the dipole operator, and *(t) is the wavefunction 
of the system obtained by solving the time-dependent Schrodinger equation (TDSE) 
with the Hamiltonian H{t): 

id<^{t)/dt = H{t)<^(t). (1.3) 

We will first focus on the single-active-electron approximation (see section 1.8 for 
the multielectron case). This approximation assumes that only one electron feels the 
laser field - the one that is liberated via strong-field ionization and subsequently re- 
collides with the parent ion. All other electrons are frozen in the ion, unaffected by 
the laser field. The Hamiltonian of our system in the single-active-electron approxi- 
mation is 

7J(t)=pV2 + C/(f) + Vi(t), (1.4) 

where p = — iVr is the momentum operator, [/(f ) describes the interaction of the 
electron with the ionic core, and Vz, (t) describes the interaction between the electron 
and the laser field. In the dipole approximation and in the length gauge, V^ (t) = 
— d • F(f) = f • F(i) (see Chapter 8 to learn about different gauges or read Section 
2.2.4 in the excellent book by |Grynberg et a/.| ( [2010[ l for a more detailed discussion). 
Formally, the solution of the Schrodinger equation ( |1.3[ l can be written in the inte- 
gral form (see e.g. [Smirnova et a/.|p007b) for a simple derivation): 



/■*..„ 
!*(*)) = -i / dt' U{t,t')VL{t')Ua{t' ,ta)\g) + Ua{t,to)\g), 

■'to 



(1.5) 



where the ket- vector \g) represents the wavefunction of the electron in the ground 
state at initial time t — Iq, U{t, t') is the full propagator, while Uo{t' ,to) is the field- 
free propagator. The propagators are the operators that describe the time evolution 
of the wavefunction. The propagator tJo{t' ,to) governs the electron dynamics from 
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time to to time t' without the laser field, and is determined by the following equations: 

idUo{t,to)/dt^HQUoit,to), (1.6) 

Uo{to,to) = l, (1.7) 

Ho = p^/2 + Uir). (1.8) 



Symbolically, the solution of Eq. \l.6\ can be written in the compact form 

Uo{t',to) = e-'^^i^°^^'>''^, (1.9) 

where the integral is time-ordered, that is, the contribution of later times to the evo- 
lution follows the contribution of the earlier times. 

The full propagator U{t, t') governs the electron dynamics from time t' to the ob- 
servation time t, driven by the combined action of the laser field and of the ionic core 
potential U{i). It is given by 

idU{t,t')/dt = HU(t,t'), (1.10) 

U{t,t') = e-'-^''"^^^'^^, (1.11) 

[/(i',t') = l. (1.12) 

The propagation without the laser field is straightforward. Denoting the ground state 
energy Eg — —Ip (ionization potential) and the stationary ground state wavef unction 
*g(r) — (r|.g), we have: 

*g(r,t') = t/o(i',to)*<,(r) = e^^-(*'-*°)*g(r). (1.13) 

The full propagator U{t, t'), on the other hand, is just as hard to find as the solution 
of the original equation \\3\ . The advantage of the integral expression Eq. \\.5\ 
is that making meaningful approximations is technically easier and physically more 
transparent. 

Remembering that the laser field is strong, we can try to neglect the ionic poten- 
tial in the full propagator. In this case the electron is free from time t' to time t. 
Its motion is only alTected by the laser field and is described by the Hamiltonian 
Hv{^) = P^/2 4- Vl(*)- The corresponding approximation is called the Strong 
Field Approximation (SEA), and the propagator corresponding to -ffy(f) is often 
called the Volkov propagator. The main advantage of the SEA is that the Volkov 
propagator can be found analytic ally. In the length gauge used here, the result of act- 
ing with the Volkov propagator Uyi't, t') on the plane wave with kinetic momentum 
k(i') = p + A(t')is 

Uv{t,t')\p + A(i')) = e-'^-(P^*'*'V + A(i)), 
(r|p + A(.)) = ^e^[P+-(*)l-r, 

Sv{p,t, t') = 1 y d^P + A(C)]'. (1.14) 
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That is, the plane wave with the kinetic momentum k(t') = p + A(t') turns into a 
plane wave with the kinetic momentum k(t) = p + A(i) and accumulates the phase 
Sv (p, t, t') on the way. 
The Eqs. ( |1.14p define the VoUcov function 

v]>^l'- f.f'\ _ 1 -»Sv(p,t,t') »[p+A(t)]-r 

P^''^" (2^)3/2^ 

Formally, the Volkov function is an eigenstate of the time-periodic Hamiltonian. It 
provides the quantum-mechanical description of the behavior of the free electron in 
the laser field. The coordinate part of the Volkov function is a plane wave, and these 
plane waves form a complete basis at each moment of time: 

i= /'dp|p + A(t))(p + A(i)l. (1.15) 

Within the SFA, Eq.([T3) takes the form 

\^{t)) = -i f dt'Uv{t,t')VL{t')Uoit',to)\g) + Uo{t,to)\g), (1.16) 

Jto 

and can be solved analytically. The first term describes ionization, the second term 
describes the evolution of the non-ionized part of the electron wavefunction. 

Thus, it is natural to associate t' with the time when ionization is initiated: before 
t' the electron is bound, after t' the electron is becoming free. Substituting Eq. (1 . 1 6 1 
into Eq. (fL2l yields: 

D(f)~-i(f/o(t,io)(7!di / dt'Uv{t,t')VL{t')Uo{t',to)\g)+c.c. (1.17) 

Jto 

Here we have assumed that there is no permanent dipole in the ground state and 
that the contribution of the continuum-continuum transitions to the dipole is negli- 
gible. The latter assumption is fine as long as ionization is weak. Thus, the dipole 
in Eq. (fTTTTl is evaluated between the bound and the continuum components of the 
same wavefunction. 

The propagator t/y (f , t') is known when it acts on the Volkov states. Thus, we in- 
troduce the identity operator resolved on the Volkov states, Eq. ifTTTS), into Eq. ( |1.17[ l: 
ft 

D(f) = -i{g\d\ / df'e'-^"'*"*) x 

Jto 

X fdpUv{t,t')\p + Ait')){p + A{t')\VL{t')\g) + c.c. . (1.18) 



Finally, remembering that VL(t) = — d ■ F(t), we re-write Eq. ( 1.18 i in the compact 
form: 

D(f)=i f dt' f dpd*{p + A{t))e-'^'-P^*-*"^F{t')d{p + A{t')) + c.c., (1.19) 

where we have introduced the dipole matrix elements d(p + A(i)) of the transitions 
between the ground state and the plane wave continuum, 

d(p + A(t)) = (p + A(i)ldl5). (1.20) 
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The phase 

1 ft 

S{p, t, t') =- I [p + A{T)fdT + Ip{t - t') 



(1.21) 



is often referred to as action, and we will use this term below, even though, strictly 
speaking, it is only the energy part of the full classical action. 

It is convenient to re- write the equation ( 1.19^ for the harmonic d ipole D(i) by 
evaluating the integral over t' by parts (see e.g. Gribakin and Kuchiev ( 1997 i, Becker 
|era/.| ( |2002b) and Appendix A): 



ft , 

/ dt'e~*^(P'*'*)F(i')d(P + A(i')) 
Jta 

rt 

= / di'e~'^(P'*'*)T(p + A(i')), 

Jto 

T(p + A(f')) = 



(P + A(t'))' 



+ Ip 



A(t')|<7), 



(1.22) 



(1.23) 



where (p + A(t)\g) is a Fourier transform of the ground state \g), T(p) reflects the 
dependence of ionization rate on the angular structure of the ground state. Equation 
(fTTT9} takes the following form: 

D(f) = i f dt' /'dpd*(p + A(f))e"''^(P'*'*')T(p + A(i')) + c.c., (1.24) 

Jto J 

The harmonic spectrum I{Nuj) can be obtained from the Fourier transform of 

D(t): 



I{Nuj) oc {Nujf\D{Niu)\'^ , 



DiNcu) = 



/ 



rffe^^"*D(f). 



(1.25) 



Note that ^(p, i, i') is large and the integrand is a highly oscillating function, which 
is an advantage for the analytical evaluation of this integral. The analytical approach 
l |Lewenstein et al.\^l994\ ) is based on the saddle point method (see Appendix B), 
which is the mathematical tool for evaluating integrals from fast-oscillating functions. 
It provides the physical picture of high harmonic generation as a three step process 
involving ionization, propagation and recombination ( jlvanov et a/.| ( |199"6| l). It also 
supplies the time-energy mapping ( |Lein| ( p005) ; [Baker et al. ( 2006) ) crucial for at- 
tosecond imaging, and it is the basis for the extension of the above approach beyond 
the SFA and beyond the single-active-electron approximation (see e.g. |Smirnova| 
|eFar] ( |2009al l). 

Let us now focus on the analytical saddle point approach to HHG. 
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1.4 

The Lewenstein model: Saddle point equations for HHG 

The goal of this section is to evaluate the integral equations l |1.24|1.25| l using the 
saddle point method (see Appendix B). We need to find saddle points for all three 
integration variables t' , t and p, i.e. points where the rapidly changing phase of the 
integrand has zero derivatives with respect to all integration variables. 

There are two ways to deal with the integrals Eqs. ( |1.24|1.25) . First, one can treat 
them as multi-dimensional integral, i.e. one finds the saddle points for all the integra- 
tion variables 'in parallel', and then one follows the multi-dimensional saddle point 
approach to deal with the whole multi-dimensional integral 'at once'. 

One can also take a different route and evaluate the multiple integrals ( |1.24|1.25| 
step by step, sequentially. First, we find the saddle points ti for the integral over t' 
from the saddle point equation: 

where the phase S is given by Eq. ( |l.2l[ . We then evaluate the integral over t' treating 
it as a one-dimensional integral, with p and t entering as fixed parameters. 

Next, we move to the integral over p. Dealing with its saddle points, we should 
keep in mind that the saddle points of the previous integral t' = ti = ti{p, t) depend 

°" p- ^ ^ 0' " = ^' y^ ^■ 



Fortunately, thanks to Eq. (1 .26 i, the explicit dependence of ti (p, t) on p does not 



affect the position of the saddle points for the p-integral: 

dS{U,Y>,t) _ dS{t,,p,t) ^ dS{U,p,t) dt, ^ dS{ti,p,t) ^^ ^^ ^7) 

dpa dpa dti dpa dpa 

Note that the integral over p is multi-dimensional, which leads to a slightly different 
form of the pre-exponential factor (prefactor) involving Hessian(see Appendix B). 

Finally, we deal with the integral over t. Here, again, the saddle points ps{ti,t) 
depend on t: ^°^'° / 0. But once again the explicit dependence of ps{ti,t) on t 
does not affect the position of the saddle points thanks to Eq. (fL27l: 

dS{t,,ps,t) ^ dS{U,ps,t) dS{ti,ps,t) dpa ^ dS{ti,ps,t) ^ ^ ^g-, 

dt dt dpa dt dt 

The fact that both routes yield the same saddle point equations is, of course, not 
surprising - one should not get different answers depending on how the integral is 
evaluated. 

Using Eq. ( |1.21[ , we obtain the explicit form of the Eqs. ( |1.26|1.27|1.28^ , which 
define the saddle points ti, ps, tr'. 

[P^iAM!+j^.O, (1.29) 

[ps + Ait' )] dt' =Q, (1.30) 

lHi±AM! + j^ = ^,. (1.31) 
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Figure 1.2 Contour of the time integration in tiie action S. Ionization occurs from the 
complex time u to the real time i^. 



Here, ps is the electron drift (canonical) momentum, ki(i) = p^ + A(f) is the kinetic 
momentum (the instantaneous electron velocity, up to the electron mass). The trajec- 
tories that satisfy the Eqs. ( |1.29|1.30|1.3ip are known as quantum orbits, see Salieresj 
|eFar] ( |200T) , |Kopold et a/."] (|2002 1, Beck er era/. | ( |2002a| l. 

Equation ( |1.30[ l requires that the electron returns to the parent ion - the pre- 
requisite for recombination. Indeed, the time integral of the electron velocity yields 
the electron displacement from ti to tr. Thus, Eq. (fL30j dictates that the displace- 
ment is equal to zero. 

Whereas Eq. (fOTJ describes energy conservation during recombination, Eq. (fL29j 
describes tunnelling. It shows that the electron's kinetic energy at t.^ is negative, its 
velocity ks(ii) = Ps -I- A(fi) is complex, and hence t^ = ^ + ^l is also complex - 
the hallmarks of the tunnelling process. 

The time ti can be identified with the moment when the electron enters the harrier, 
see Fig. |1.2| Its real part will then correspond to the time when the electron exits the 
harrier. The origin of this concept will be explained in the next section. 

The electron displacement during this ' under- the-barrier' motion from ti to Re(ii) 
is, in general, complex. Whether we like it or not, it yields a complex coordinate of 
'exit' Vex ~ v'ex + iv'lx at Re(f j) = t'i (see e.g. Torlina and Smirnova 1 2012 ): 



/: 



[p + A(t')]dt' = r'ex + ire 



(1.32) 



As a result, Eqs. ( |1.30[ l and ( |1.31[ cannot be satisfied unless p or tr are complex. In- 
deed, tr must be complex to compensate for the imaginary displacement accumulated 
under the barrier. However, the energy conservation condition in Eq. ([T3T) dictates 
that ks(fr-) = p^ + A{tr) is real at the moment of recombination. Therefore p., must 
also be complex to compensate for the imaginary part of A(tr). 
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Thus, we are forced to conclude that, in contrast to the classical trajectories of the 
simple man model, the quantum orbits are the trajectories with complex canonical 
momenta, complex velocities, and complex displacements. These trajectories evolve 
in complex time. The only quantity that is required to be real is the one we measure 
- the energy of the emitted photon, see Eq. (fOTl. Later in this chapter, we will see 
when and how one can replace these trajectories with a different set of trajectories 
that do not involve complex canonical momenta and therefore better correspond to 
the classical picture. But for the moment, let us deal with the problem at hand. 

For a linearly polarized field, it is convenient to rewrite the Eqs. ( |1.29|1.30|1.31[ l in 
terms of electron momenta parallel, Ps \\, and perpendicular, p^ j^, to the polarization 



vector of the laser field: 

l2 



\Ps,\\+Mt^)] 



+ Ip,cff = 0, 



i: 



[p,«+A{t')]dt' = 0, 



dt' = 0, 



\PsA\+Mtr)f 



+ h 



,Gff 



= iVLJ, 



(1.33) 
(1.34) 

(1.35) 



where we have introduced an "elfective" ionization potential: /p. off = Ip + p^ j_/2. 
Equations (|1.34[l dictate that the stationary perpendicular canonical momentum is 



equal to zero for the linearly polarized field, pg± 
Eqs. ( |1.33|1.34|1.33| l reduce to: 

ftr 

I b,_||+A(i')]dt' = 0, 



\PsA\+Mtr)Y 



+ Ip = Nuj. 



and hence Jp^cff = Ip ■ Then, 

(1.36) 
(1.37) 

(1.38) 



Separating the real and the imaginary parts in Eqs. ( |1.36|[L37]|1.38^ , we obtain six 
equations for six unknowns: ti — t[+ it", U = tr+ it", p^ u = p' + ip" . Our goal is 
to solve these equations for each harmonic order A'^. Here is one way to do it, which 
we find simple and visually appealing. 

First, we use the Eqs. ( 1.36|l.38 i to express all variables via the real, tj,, and the 
imaginary, i", return times. This can be done analytically. Second, we substitute the 
result into the real part and the imaginary part of Eq. ( |1.37^ : 



Fi{N,t'r,tr) 
F2{N,t'r,tr) 



Re 



Im 



[Ps 



A{t')]dt' 
j\ps.,\\+A{t')]dt' 



(1.39) 



(1.40) 



Third, we solve the Eqs. ( |1.39|1.40| l to find the only two remaining unknowns: the 
real, fj,, and the imaginary, t", return times. While the Eqs. (1.39 1.40 1 cannot be 
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solved analytically, dealing with two equations is much easier than dealing with the 
original six. 

Solving the Eqs. ( |1.39|1.40} means that we need to find the minima of the two- 
dimensional surface F{N, t'j.,tr), defined in the plane of the real, t'r, and the imagi- 
nary, f", return times: 

F{N,t'r,t';) = [Fi{N,t'r,tr)? + [F2{N , t'^ , t'r)? = 0. (1.41) 

These minima can be easily found numerically using the gradient method. The advan- 
tage of using -F'(A'', t'^., i") is the ability to visualize the solutions: by simply plotting 
the surface given by Eq. (fL4TJ, see Fig. |1.3[ one can examine the positions of the 
minima versus the harmonic number A'^. 

If we restrict our analysis to those solutions that lie within the same cycle of the 
laser field as the moment of ionization, Re(ij), we will find two stationary solutions 
for each harmonic number TV. These solutions are discussed in detail in the next 
section. They correspond to two families of quantum orbits, called the 'short' and 
the 'long' trajectories. The trajectories merge for the largest possible return energies, 
i.e. near the cut-ofF of the harmonic spectrum. 

There are also solutions that lie outside the laser cycle during which the electron 
was 'born' into the continuum. These 'super-long' trajectories describe second, third, 
and higher-order returns of the electron to the origin. In typical experimental con- 
ditions, their contribution to the high harmonic emission is negligible thanks to the 
macroscopic effects - very long trajectories do not phase match well (see e.g. Salieresj 



|ef a/.|([200T)). Only very recently, the beautiful experiments of|Zair et al. |([2008 1 have 



been able to clearly resolve the contribution of these trajectories, and even identify 
their interference with the contribution from the long and the short trajectories. 

The stationary phase method for the integral over the return time t breaks down 
when these two stationary points merge and the second derivative of the action with 
respect to the return time is equal to zero, d^S/dt^ — 0. At this point, one needs to 
replace the standard saddle point method with the regularization procedure, discussed 
in Appendix C. 

Outside the cut-off region, and up to a global phase factor, the saddle point method 
yields the following expression for the harmonic dipole (|1.24|1.25): 



4A/ 



D(iV^) = Yl 



2-K 

W7. 



2-R 



(27r)3/2 



Vdet(iS^',,pJ 

xd (ps + A(4-"))e ^*^"' >■ '* ^ T(ps + A(ij))e -- , (1.42) 

where the Hessian det(iS'p^,p^) appears due to the multi-dimensional nature of the 
integral over p. The sum runs over all stationary points j for M periods of the laser 
light, and the corresponding ionization and recombination times are labelled with the 
superscript j. Since there are two trajectories for each half-cycle of the laser field, 
i.e. for each ionization 'burst', and since there are 2M ionization bursts for M laser 
cycles, the number of stationary points is AM. 
The length gauge SFA presents a good approximation for short range potentials 
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[Frolov et al. | ( |2003y However, it misses polarization, the Stark shift and the depletion 
of the bound state, which can be introduced into Eq. ( |1.42^ if necessary. 

Note that the expression (fL42l can not be directly ported to long range potentials. 
Indeed, in the long range Coulomb potential the ground state has different radial struc- 
ture (compare Eqs.i 1.551 and ( 1.56 1 below) and T(ps + A(fi)) is singular exactly 
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at the saddle point [ps + A{ti)f/2 + Ip — Q. As a consequence, the saddle point 
calculations should be modified, (see e.g. |Keldysh| ( |1965) ; [Gribakin and Kuchiev| 
( |1997) ) to accommodate for the presence of such singularity. Once the singularity 
is treated correctly (Keldysh|( [1965^ ; pribakin and Kuchiev| ( |199 7')), the result is still 
incomplete and unsatisfactory, because the long range potential also affects the struc- 
ture of the continuum states, which can not be accurately represented by the Volkov 
states. Consistent treatment of long-range effects including modifications of both 
bound and continuum states can be found in |Torlina and Smirnova| ( |2bl2) ; |Kaushal| 
[and Smirnova| ( pOT3) . A practical recipe for incorporating the effects of the Coulomb 
potential is discussed in the next section and in the Outlook section. 



1.5 

Analysis of the complex trajectories 

Let us now show how the above method of finding the saddle points works for a 
linearly polarized laser field F — Fq cos(a;i), which corresponds to the vector po- 
tential A — — Aq sin(ajt). We shall introduce the dimensionless variables pi = 

Re(Ps_||)/^0, P2 = Im(Ps,Jl)Mo, (t>i = ^ti = 01 + ■i(t>i, 4>r = Lotr = (t>'r + i4>r , 

^2 = /p/(2C/p), 'yl = [Nuj - Ip)/ {2Up) 
In terms of these variables, Eqs. l |1.39|1.40) for the linearly polarized field yield: 

Fl = pi(0r - 0j) -P2(0r - 0") -COs((^i)cOsh(</>f) 

+ cosh{(l>r) cos{<t>r) = 0, (1.43) 

F2 = pi{(l>r - (l>i)+P2{<t>r ~<t>i) +sin{(f)',)smh{(f)i) 

- sinh((^") sm{(j>'r) = 0. (1.44) 

The real and the imaginary parts of Eq. ([L38j allow us to express the real, pi, and the 
imaginary, p2, components of the canonical momentum via the real and the imaginary 
parts of the return time (for above threshold harmonics): 

pi = COsh{(t>r) sm{(j)r) + JN, (1-45) 

P2 = sinh{<l)r) cos{(j>r) ■ (1-46) 

The real and the imaginary parts of Eq. ( |1.36^ , 

pi = cosh(<^") sin(0j), (1.47) 

P2 + 7 = smh{(f)'^) cos{(j>i), (1.48) 

allow us to express the real, t/i^, and the imaginary, 0", ionization times via pi and 
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Figure 1.3 Surface in Eq 



Imaginary time of return, units of laser cycle 

for /„ 



1 
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1 



Imaginary time of return, units of laser cyole 



1.41 



15.6 eV, 7 = 1.3 ■ 10" W/cm^, Hoj = 1.5 eV. Left 
panel: N = 11; two minima corresponding to short (t/ij, ~ 0.4) and to long (<^'^ ~ 1) 
trajectories, respectively. Right panel: N = 27; the two minima corresponding to short and 
to long trajectories are merging together. 



P2'- 



0^ = arcsin(V(P-L>)/2), 
0" = arcosh(A/(P + D)/2), 



(1.49) 
(1.50) 



where 



P^pI+1^ + 1, 



7 = 7 + P2- 
Now we can use our recipe: 



(1.51) 



using 



/)'/intoEqs. (1.43 1.44 1; 



Fi + f| in the plane of the real and the imaginary return 



• Pick a grid of values (f>r, 4>r in the complex plane of the return times . 

• Pick a point (f>r, 4>r and calculate pi, p2 using Eqs. ( 1.45 1.46 1 and <; 
Eqs. ( |1.49|1.50| i; 

• Substitute pi, p2, (j)' 

• Plot the function F 
times; 

• Look for the minima, see Fig. |1.3| 

Instead of reading out the solutions from the graph, one can find the minima using 
the gradient method. An alternative algorithm using the same ideas is described in 
Appendix D. 

The imaginary and the real return times (Fig. |1.4| right panel) define the integration 
contour in the complex plane: only along this contour the energy of return and there- 
fore the energy of the emitted photon are real. This energy is shown in Fig. |1.4| (left 
panel) vs the real component of the return time for typical experimental conditions. 

The cut-off (maximal energy) corresponds to about 3.17 Up + 1.32 Ip, see lLeweiT] 
[stein erar| ( | 1994) . There are two different trajectories returning at different times that 
lead to the same re-collision energy. Those returning earlier correspond to shorter 
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Figure 1.4 Left panel: Emission energy, E{tr) 



0.2 0.4 0.6 0.8 

Real time of return, units of laser cvcle 

= Afw, vs real time of return for 



Ip = 15.6 eV, 7 = 1.3 ■ 10^'' W/cm^. Right panel: Imaginary time of return vs real time of 
return. The solution diverges in the cut-off region. The thick blue line schematically shows 
the desired outcome of the regularization procedure. 



excursion and are called 'short trajectories', those returning later are called 'long tra- 
jectories' as they correspond to larger excursions and longer travel times. 

Fortunately for attosecond imaging, the contributions of the long and the short 
trajectories to the harmonic emission separate in the macroscopic response: the har- 
monic light diverges differently for those trajectories, and thus the signals coming 
from short and long trajectories can be collected separately. As a result, each har- 
monic TV can be associated with a particular time delay between ionization and re- 
combination, t'r—t[, and therefore each harmonic takes a snapshot of the recombining 
system at a particular moment of time. This time-energy mapping, see ( |Lein| ( (2005^ , 
[Baker et a/.| ( |2006") , |Shafir et a/.|P012^ ), is the basis for attosecond time resolution 
in high harmonic spectroscopy. 

As mentioned in the previous section, the stationary phase (saddle point) method 
for the integral over return times t breaks down near the cut-off, where the two station- 
ary points (short and long trajectories) begin to coalesce and the second derivative 
of the phase S with respect to the return time is equal to zero, d^S/dt^ = 0. The 
regularization of the solutions in the cut-off region is discussed in Appendix C. Here, 
we shall proceed with the analysis of the stationary phase equations and turn to the 
ionization times. 

The concept of ionization time together with the semiclassical (trajectory) perspec- 
tive on ionization has been first introduced by V. Popov and co-workers (see |Pereio^ 
|mov et al. | ( [19661 [T967) ; rPerelomov and Popov| ( fT967) ; |Popov et al. | ( fT968) ). Just like 
in the Lewenstein model described above, the concept of trajectories arises from the 
application of the saddle point method to the integral describing ionization rj 

3) Eq. |l.52| corresponds to the length-gauge SFA result for ionization. Eq. jl.52| also results from 
the PPT approach under the approximation, substituting the laser-dressed hound wave-function by the 
field-free ground state. Thus, the PPT approach allows one to identify the approximations in Eq. jl.52^ 
for short-range potentials. The SFA is inaccurate, because the Vokov states are not sufficiently accurate 
even for short-range potentials. 
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./o 



(1.52) 



Here the upper limit of the integral in the action ^(p, T, t') Eq.( |l.2l| is the real time 
T, at which the liberated photoelectron is observed (detected), function T(p + A(i')) 
(see |1.23^ contains the Fourier transform of the bound state wave-function. The sad- 
dle point method applied to Eq.( |1.52p yields the ionization amplitude: 

nl/2 



.(P,i) 



2-K 



'^^P'^'''^T(p + A(iO), 



where ij is the complex saddle point given by the condition 

ip±4M! + 7,.o. 



(1.53) 



(1.54) 



The saddle point method selects specific moments of time f j when ionization occurs. 
At these times the instantaneous electron energy 2 ^^ equal to the energy of 

the ground state and therefore the instantaneous momentum of the electron hits the 
pole of the bound state wave-function in the momentum space. In the vicinity of the 
pole, the wave-function in the momentum space is determined by the asymptotic part 
of the wave-function in the coordinate space. Thus, in contrast to one-photon ioniza- 
tion, which probes the bound wave function near the core, the strong field ionization 
probes the asymptotic part of the bound wave function: 



(r|p) ^ C„,^3/2 e_ 



-Yu 



r 



(1.55) 



Equation (fT35l restricts our analysis to short range potentials (see [Torlina and| 

[Smirnova] pO 1 2[ l for consistent analytical treatment of strong-field ionization from a 

long range (Coulomb) potential), k = 

angular str ucture of the bound state. For Coulomb potential 

expression (1.55 1 must be multiplied by (nr)'^''^: 



2Ip, C^i is a constant, Yi^{^) reflects the 



-Q/r, the asymptotic 



{Ag) - Cf^in 



3/2 



■(«r)0'-y, 



.(;)■ 



(1.56) 



Evaluating the Fourier transform of (1.55 1 we obtain explicit expression for T(p) 
(see [Perelomov et al. \ (JT966J) : 



T(P) 



I CulYlmi — ) 

■K I p 



(1.57) 



Perelo- 



Evaluation of the spherical function y/m( p) ^t the pole p — iJAc yields (see 

|mov e?ar| ( |1966[ l and also |Barth and Smirnova|p01 1[|2013) for circularly polarized 

fields) : 



Mm(^)lp=±iK — Qm ( I 



•p' 



imrj>p 



Clr^ 



{21 + 1){1 + \rn\ 



2l™ljmiV 47r(;-|m|)! 



(1.58) 
(1.59) 
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where p± is the electron transverse momentum at the detector and (f>p is the azimuthal 
angle of the electron momentum at the detector. Taking into account that S'/. f. = 
(p + A{ti))A'{ti) can be written as 



II = -sign(F(tj)), 



(1.60) 
(1.61) 



we obtain the final expression for the amplitude of single ionization burst at complex 
time ti{p), specified by the final momentum of the electron p: 



a-ionip) = 2C 



fiF{t, 



P_L 



'g-iS(p,T,t.(p))gjm0p 



(1.62) 
(1.63) 



Note that both the real and the imaginary component of ti (p) depend on laser param- 
eters. Therefore, the dependence of the pre-exponential factor (prefactor) in aio„(p) 
on the laser field is not simply F^^''^ as illustrated below for the ionization at the 
maximum of the laser cycle. Indeed, the sub-cycle dynamics in the prefactor is much 
slower than in the exponent, thus one can also use a simpler expression for 5". ^. cor- 
responding to its value at the maximum of the laser field ( [Perelomov et a/.| ( |l966^ ): 

<,.=-^^^^. (1.64) 



Omitting the sub-cycle dynamics in the prefactor we obtain a simpler expression for 
the amplitude of single ionization burst at the time ti(p), consistent with the one 
derived by [Perelomov et a/T| ( |1966) : 

nl/2 



(p) = 2C 



-^LO 



K\/l + 7^ 



m' 



„-JS(p,T,t,(p))-(-i 



(1.65) 



At the same level of approximation, i.e. neglecting the sub-cycle dynamics in the 
prefactor, the effects of the Coulomb potential are incorporated by simply adding the 

/„ a\Q/K 

factor ^ 



aioTi(p) =2C 






-7CJ 



=\/n^ 






-iS{p,T,ti{ji))+iva(f>p 



(1. 



The sub-cycle Coulomb effects are derived in |Torlina and Smirnova|p012^ . Note 
that in the rigorous analysis within the analytical R-matrix (ARM) approach, which 
consistently treats the Coulomb effects both in bound and continuum states ( [Torlina| 
and Smirnova (2012 1; Kaushal and Smirnova (2013 ), the pole in T(p -I- A(ti)) does 
not appear, because the radial integration is removed due to the use of the Bloch 
operator ( [Bloch|p957^ ). Therefore, it removes all technical aspects and additional 
terms associated with the presence and the strength of the pole. 

Note that the expression for the induced dipole, Eq. (fL42l, contains terms that look 
very much like the ionization amplitude Eq. (fT33}. This observation is important. 
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Figure 1.5 Left panel: Real and imaginary ionization times vs real return time for 
7p = 15.6 eV, 7 = 1.3 ■ lO^"* W/cm^, huj = 1.5 eV. Right panel: Cartoon illustrating the 
ionization window. Ionization occurs around the field maximum within an approximately 
250 as time window (corresponding to the maximum value of the real ionization time). 



as it suggests the connection of the harmonic response to ionization, as in the simple 
man model. However, the story is more subtle: the stationary momenta ps in the 
harmonic dipole are complex-valued, while here they are real observable quantities. 

The integral Eq.( |1.52[ l has been extensively studied by Keldysh, Popov, Perelomov, 
Terent'ev, and many others. The semiclassical picture in |Perelomov et ar| ( |1966| 
1 1967) ; [Perelomov and Popov] ( |1967[ l; [Popov et a/.| ( |1968) , enabled by the applica- 
tion of saddle point method, shows that strong-field ionization can be understood as 
tunnelling through the oscillating barrier created by the laser field. The tunneUng 
picture clarifies the sensitivity of strong field ionization to the asymptotic 'tail' of the 
bound wave-function (see Eq.( |1.55|1.56) ), since it is this asymptotic part that 'leaks' 
through the barrier. The modulus of the ionization amplitude is associated with the 
imaginary part of the action S in Eq. ( |1.53|[L62]|1.65[[l.66) . This imaginary part 
is only accumulated from ti to t', since in the photoionization problem the canon- 
ical momentum registered at the detector is real and the integration over time also 
proceeds along the real time axis between t[ and the observation time t. 

This is why the complex saddle point ti is associated with the time at which the 
electron enters the classically forbidden region - the tunnelling barrier - while the real 
part of the complex saddle point t[, after which changes to the ionization amplitude 
stofTj is associated with the time of exit from the classically forbidden 'under-the- 
barrier' region. The same reasoning can be extended to the ionization times arising 
within the semiclassical picture of harmonic generation, see Fig. lfT3J. However, 
the ionization times in high harmonic generation are somewhat different due to the 
fact that -ps are complex-valued. In the next section we will consider the connection 
between these two times. 

The imaginary ionization time defines the ionization probability. Since the imagi- 



4) Rigorously, this statement is only true for short range potentials. Long-range electron-core interactions 
lead to additional modifications of the ionization amplitude after t'^ jTorlina and SmirnovaH2Q12| , |Tor-| 
|linaefa].|l2013i). 
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Figure 1.6 Left panel: Real canonical momentum vs real return time for Ip = 15.6 eV, 
/ = 1.3 • 10^* W/cm^, huj = 1.5 eV. Right panel: Imaginary canonical momentum vs real 
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nary component of the ionization time is larger for short trajectories, these trajectories 
have a lower chance of being launched compared to the long ones. The range over 
which the real part of the ionization time changes within the quarter-cycle defines 
the duration of the 'ionization window'. Typically, for a A ~ 800 nm driving laser 
field and a laser intensity of / ~ 10^ W/cm^, the ionization times (their real part) are 
spread within ~250 attoseconds around the instantaneous maximum of the laser field 
(see Fig.[T3J. Thus, strong-field ionization is an intrinsic attosecond process. Note 
that the quantum 'ionization window' is shorter than the classical one (see Fig.fTTl, 
as according to the classical simple man picture ionization happens at any phase of 
the laser field. 

Figure [L6] shows the saddle point solutions for the electron canonical momentum. 

In photoionization, the electron canonical momentum is always real, since it is 
the observable registered at the detector. In contrast, in harmonic generation the ob- 
servable registered at the detector is the emitted photon, and hence it is the photon 
energy that must be real. As a result, the electron canonical momenta in HHG are 
complex. Electrons on long trajectories have a very small imaginary canonical mo- 
mentum. Therefore, it is a very good approximation to associate long trajectories 
with photoelectrons. Note that the maximum of the real canonical momentum is 
about pmax — Aq. In the photoelectron perspective pmax corresponds to an energy 
of 2Up at the detector - the cut-off energy for the so-called direct photoelectrons, i.e. 
those that have not substantially changed their momentum after ionization. 

The imaginary part of the canonical momentum can be quite large for short tra- 
jectories. The complex-valued solutions, not only for the ionization times, but also 
for the recombination times and the electron canonical momenta, challenge our un- 
derstanding of the underlying physical picture of harmonic generation. If the first 
step of high harmonic generation is ionization, then why do these liberated electrons 
have complex canonical momenta? Does this mean that these electrons have not been 
ionized? Can we factorize the harmonic dipole into ionization, propagation and re- 
combination? The next section explores this opportunity. 
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1.6 

Factorization of the HHG dipole: simple man on a complex plane 

Having derived the analytical expressions for the HHG dipole, can we identify the 
simple man model in it, within the consistent quantum approach? To do this, we 
need to factorize the harmonic dipole into the three steps: ionization, propagation and 
recombination. That is, we have to re-write the dipole as a product of the ionization 
amplitude, the propagation amplitude and the recombination amplitude. 

Such factorization of the harmonic dipole is not just curiosity driven. It is important 
for extending the modelling of harmonic emission to complex systems. Once the 
three steps are identified, the respective amplitudes can be imported from different 
approaches, tailored to calculate specifically ionization or recombination in complex 
systems. 

The factorization of the harmonic dipole runs into two types of problems: technical 
and conceptual. The technical problems arise from the fact that the original three- 
step (simple man) model is formulated in the time domain. The three processes - 
ionization, propagation and recombination - are the sequence of subsequent time- 
correlated events. The harmonic spectrum formally corresponds to the harmonic 
dipole in the frequency domain, where the three processes become entangled: recall 
the contribution of different quantum trajectories to the same photon energy. Thus, 
rigorous factorization in the frequency domain is only possible in the cut-off region, 
where short and long solutions merge, see Frolov efa/.|(|2009,lp*| 



The conceptual problem is due to the complex canonical momentum of the elec- 
tron responsible for HHG. Ionization in terms of creating photoelectrons with real 
canonical momenta does not appear to fit into the HHG picture. Can we build an al- 
ternative model of HHG based entirely on photoelectrons, i.e. those electrons which 
are indeed ionized at the first step? 

Let us address these issues step by step, starting with the factorization of the har- 
monic dipole in the frequency domain ([Frolov et a/.|(|2009 1; Kuchiev and Ostrovskyj 



| [T999) ; [Morishita et a/.| pOOSj i) and the time domain ( |Ivanov et dl] (1996} ). The 



former involves the factorization of Eq. ( |1.25) , the latter factorizes Eq. ( |1.24[ l. 

1.6.1 

Factorization of the HHG dipole in the frequency domain 

To re-write the harmonic dipole in the semi-factorized form, we can take Eq. (fL42l 
and split the action integral S that enters the phase of this expression into the follow- 
ing three time intervals: from ti to t'^, from ^ to ij,, and from t'^ to tr (see Fig. |l.2[ . 
Then we can identify the group of terms that looks like the ionization amplitude sim- 

5) Note that the quantitative rescattering theory (see |Le et al. |J2QQ9| ) postulates that one can factor out the 
recombination step in the frequency domain harmonic dipole. This postulate is supported by the results 
of numerical simulations demonstrating approximate factorization in the cut-off region, see jMorishitaj 
[rt^|20O8). 
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ilar to that given by Eq. (|T33j, 



O-ionyPs, ti) — 



2-R 
W7, 



1/2 



-iS{ps,t[,ti) 



T{ps + A{U)). 



(1.67) 



The ionization amplitude is associated with the first time interval, from ti to i^, and 
only the part of the action integral from ti to the real time axis t'^ = Re(ti) enters 
this amplitude. For a short-range potential (neglecting the sub-cycle effects in the 
prefactor): 



aion(Ps) 



2C 



-JUJ 



Kv/TT^ 



Ps± 



-iS{ps,T,ti)+im4 



(1. 



Constant C is specified in Eq. ( |1.63[ l. The momentum ps is given by the full set of 
saddle point conditions for ti,tr, and pst 



[Ps + Hk)? 



Ip^O, 



[ps+A{t')]dt' = 0, 



[ps+A{tr) 



+ Ip 



Nuj. 



(1.69) 



Note that ti in HHG and and ij in ionization are different, that is why in Eq. |1.65| we 
use 7 and p, whereas in Eq. ( |1.70[ l we use 7 (see Eq |1.51[ l and ps. If imaginary part 
of ps is equal to zero, then 7 = 7. For Coulomb potential (neglecting the sub-cycle 
effects in the prefactor): 



(p.) = 2C 






-70; 



^VTP 



(PsAV^ ^-iS(p,,T,U)+ini<t,^„ J.J .^Q^ 



Now consider the next time interval, from i'^ to t'^. The prefactor arising from saddle 
point integration over the electron momenta p leads to the term 

(27r)3/2 _ (27r)3/2 



Vdet(zS^',,pJ (i(i,-t,))3/2- 



(1.71) 



This term describes the free spreading of the electron wavepacket between ti and tr 
Thus, we associate the following group of terms with the propagation amplitude: 



CLprop \Ps ^tr jti 



(27r 



,3/2 



(i(f.-i,))3/2 



-iS(ps,t'r,t-) 



(1.72) 



Note that the denominator includes the complex-valued times f j and tr- 

Finally, the recombination amplitude is represented by the recombination matrix 
element d* [ps + A(tr)) and can be associated with the following group of terms: 



= (Ps,ti 



2n 



^s'U 



1/2 



-iS{ps,tr-,t'^)+iNutr. ■,* 



d*{jps+Htr)), (1.73) 
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where S't^tr ~ ^\/'2{Nuj — Ip)F(tr) for a linearly polarized field. As a result, the 
total dipole is formally written as 

4A/ 

where the index j labels the saddle points. However, in contrast to photoelectrons, 
the electrons involved in HHG have complex canonical momenta ps. Therefore, the 
imaginary part of the action is accumulated not only 'under the barrier', from ti to t'^, 
but also all the way between t[ and tr- Thus, factoring out ionization as the first step 
of HHG is not that convincing. Similarly, the recombination step involves not only 
the recombination dipole, but also the possible change in the amplitude due to the 
imaginary contribution to the action between fj. and tr- Thus, while we can formally 
associate several groups of terms in the harmonic dipole (fL74l with amplitudes of 
ionization, propagation, and recombination, the complex-valued electron momenta 
make such identification somewhat stretched. 

An additional point to note is that the three amplitudes are also entangled due to 
the sum over the different saddle points in Eq. ( |1.74[ l. Even if we only consider con- 
tributions of the two most important trajectories, short and long, the sum entangles 
their contributions and also mixes up the contributions from different half-cycles. 
Importantly, a finite pulse duration leads to a different mapping between the given 
harmonic number and the ionization-recombination times for each half-cycle. 

A practical approach to fac torization realized in the so-called quantitative rescat- 



tering theory (Le et al. (2009 i) is to assume that airec{'Ps{N),tr {N)) — a.rec(Nu>) 



for all j yielding: 

4M 
Ii{Nij) = a.rec{Nuj)Y,'^prop{ps,ti^\t^f^)a,oniPs,t[^''), (1.75) 

J = l 

This approximation breaks down in the following cases: 

1) In two-color orthogonally polarized fields |Morales e?ar| ( [2012^ . In this case more 
than two trajectories returning at different angles can map into the same return 
energy [Morales et a/.|p012| l. Such trajectories must correspond to different re- 
combination dipoles for different angles, violating |1. 75 1 

2) In the vicinity of the structural minimum of the recombination matrix element, or 
when the phase of the matrix element changes rapidly ( [Smirnova et al.\ ( |2009b^ , 
[Patchkovskii et a/.|pOT2l l). 

3) When the sub-cycle dynamics associated with the electron interaction with the core 
potential can not be neglected. 

These technical problems can be remedied by looking at the dipole in the time do- 
main. 
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There are several advantages of using the time-domain dipole. For starters, if we do 
not perform the Fourier transform analytically, the time tr no longer has to be com- 
plex. With the Fourier integral performed using a standard FFT routine, we can keep 
tr on the real time axis. The number of saddle-point conditions is also conveniently 
reduced to two (one of them, for the momentum ps, is in general 3-dimensional) 



[p. + Mu 



+ Ip^O, 



i: 



[ps + A(i')] dt' = 0, 



(1.76) 



with tr being the parameter, instead of the harmonic number A''. 
In the time domain, it is natural to sort the contributions to the induced dipole 

according to the corresponding ionization bursts. Then, for each half-cycle j, there 

(i) 
is a single ionization burst j at time t^' that contributes to the induced dipole as 

a function of the real return time tr, see left panel of Fig. |1.5| After saddle-point 

integration, this contribution is: 



Tt^'Hu) 



2n 



iS", 



xd*(ps 



1/2 



(27r 



,3/2 



Vdet(J5'^'^,pJ 



A(t.))e-*«(P'*'-*"')T(p.+A(ip))), 



S{ps,tr,ti)= - / [ps+A{r)YdT + Ip{tr-t,), 



(1.77) 



i3/2 



1.71 I. Just as in the frequency 



with det(iSp^^pJ = \i(tr - t\^^)\ (see also Eq. 

domain, up to a global phase factor the dipole can be written as a product of three 

amplitudes: 



,0) 



(i)N 



+0■)^ 



JJ [tr) — ^rec\'Psfir)0>prop['Ps-,tr,t^ )0'ion\Psiij^ ) 



(1.78) 



The ionization and the propagation amplitudes entering this expression are given by 
Eqs. ( |1.67|1.72) . The recombination amplitude is simply equal to the recombination 
matrix element d*(ps + A(tr)), as we have not performed the Fourier transform 
yet. Equation (fL78j is the natural mathematical formulation of the three step model, 
which is intrinsically sub-cycle. 

If we ignore multiple returns and very long trajectories, then for each tr there is 
only one ionization burst to deal with. As opposed to the frequency domain, the con- 
tributions of the long and the short trajectories from this ionization burst are not yet 
mixed - they are separated in time. This is very convenient if you need to look at the 
contribution of only the short, or only the long trajectories: it is straightforward to 



Thomas Schultz and Marc Vrakking: Attosecond and Free electron Laser Science 
Chap. 1 — 2013/4/10 — 2:48 — page 24 



24 



add a time-domain filter that would filter out the unwanted contributions. Essentially, 
this would correspond to making a window Fourier transform of the time-domain har- 
monic dipole. The inclusion of the contribution of multiple returns is rarely required 
for typical experimental conditions. 

To model the full D(ir-) one needs to model ionization, recombination and prop- 
agation separately for each half-cycle, and then collect the contributions from each 
half-cycle (each ionization burst): 

D(iO = 5^D(-'')(iO. (1.79) 



To obtain the harmonic spectrum, we have to perform the Fourier transform, which 
is convenient to do numerically using a FFT routine. There are two possible ap- 
proaches to implement the Fourier transform. 

Integration along Lewenstein's contour. In this approach, the Fourier transform 
is performed along the time contour in the complex plane tr = ij- -I- it". In this case 
the argument of the recombination dipole ps -f A(ir) remains real and so does the re- 
collision energy Erec{tr)- Since it is difficult to numerically perform an integration 
along a complex contour 



DiNiu) 



dt 



r e 



-Nu)t'' 



T>(tr 



iNbjt' 



one can use variable substitution and integrate over the real return times t'j. : 



D{Niu) = / dtr 



dtr -\E^ 
e 

dt'r 



o(tr)+IpK 



Wr. 



iNut' 



Nlo = EreciU) + Ip, 
Erecitr) = [ps + A{tr)f /2. 



(1.80) 



(1.81) 

(1.82) 
(1.83) 



The derivative i n the s quare bracket is associated with the variable substitution. 

Note that Eq. ( 1.81 i contains one approximation: the term e"' '*>- is modified 
according to the energy conservation Nlj — Erecitr) + Ip- However, the integration 
of Eq. (fLSTJ is not very convenient due to the additional effort associated with the 



dt^ 



1.4 



need to avoid the divergence of ^ in the cut-off region (see Fig. 

HHG dipole on the real time axis. To keep things simple, one can keep the half- 
cycle harmonic dipole on the real time axis: 



Yi'-^\t) = SLrec{Ps,t)aprop{p. 
Ai) 






liPsyt 



U)-, 



(1.84) 



where the saddle points p^ and tf' are given by the Eqs. (1 .76 i and the index j labels 



different solutions corresponding to the same return time t. 

In this 'real-time-axis' approach, the return time f is a parameter: we have to find 
ti and Ps for each t. This can be done using a procedure similar to that described 
in the previous section, only simpler. Specifically, we introduce the dimensionless 
variables (j) — ut and p/{F/ui) = pi + ip2- For a linearly polarized field p^ _l = 0. 
For each real we use the Eqs. ( 1.43|l.44 with (^" = and (pr = (t>: 



El (0) = Pi (0 - flii) + P20" - ' 



s(0-) cosh(0-') + cos(<;/)) = 0, (1.85) 

dn(<j!>-)sinh(<;/>-') = 0. (1.86) 
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We can now use the Eqs. ( 1.47 1.48 1 to express pi and p2 in terms of 0j and 0^ . 
Then, we build the surface F{(t>) = F^ + Ff for each (j). Next, we find the minima on 
this surface. Alternatively, we can use pi and p2 as our variables, expressing (f)[ and 
0" via pi,P2, then the minima on the surface F^ + f| will yield the real, pi, and the 
imaginary, p2, components of the canonical momentum, and then the Eqs. ( |1.49|1.50l l 
yield the corresponding ionization times. 

In this approach the divergence at the cut-off is avoided, since the divergence oc- 
curs in the complex plane of the return times when calculating the Fourier transform 
analytically using the saddle point method. The price to pay is that the recombination 
dipole has to be taken at the complex arguments ps + A(i) and the re-collision energy 
Erec (t) has an imaginary part. In practice, one can use the real part of the re-collision 
energy as the argument of the recombination dipole. If one wants to avoid this ap- 
proximation, one has to extend the recombination dipoles into the complex plane of 
the electron momenta. 

Thus, one can formally factorize the harmonic dipole in the time domain, over- 
coming the technical problems associated with the factorization. However, one has 
to keep in mind that the ionization amplitude has to be modified to include com- 
plex canonical momenta and slightly different ionizaton times. Fortunately, it does 
not lead to changes in angular factors, because T„(ps + A(ti)) remains the same. 
Indeed, both ps and ti are different in case of HHG and ionization, but the term 



Ps -I- A.(ti) = i^K, (see Eq. 1.61 1 is the same in both cases. The changes appear in 
the phase S{ps,tr,ti) and the sub-cycle core effects, i.e. everywhere where ps and 
ti contribute separately. 

The conceptual problem associated with understanding the physical meaning of 
the complex electron momenta, especially in the context of the "ionization step", 
still remains. The next section shows how, and to what extent, this problem can be 
circumvented. It introduces the photoelectron model of HHG, where the electron 
canonical momentum is restricted to the real axis. 



1.7 

The photoelectron model of HHG: the improved simple man 

In the standard simple man model, the electron motion between ionization and re- 
combination is modelled using classical trajectories. Naturally, the electron velocity, 
the ionization time, and the recombination time are all real- valued quantities. In the 
quantum description, the rigorous approach based on the saddle point method leads 
to trajectories with complex-valued momenta and complex-valued ionization and re- 
combination times. The presence of complex canonical momenta makes it difficult 
to identify the ionization step. 

The complex-valued canonical momenta and recombination times arise from the 
requirement that the electron returns exactly to its original position. Since the tun- 
nelling electron accumulates an imaginary displacement during its motion in the clas- 
sically forbidden region, the complex-valued momenta and return times must com- 
pensate for this displacement. 



25 
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■ Imaginary time 



<t Tunnel entrance ti=ti'+ti" 
t 



V 
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_L_ 



Real time 
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Time of birth units of laser t^vcle 



Tunnel exit 



Figure 1.7 Left panel: Complex ionization time for pfiotoelectrons f^.p^ vs time of birth for 
Ip = 15.6 eV, / = 1.3 ■ IQi"' W/cm2, hai = 1.5 eV. Right panel: Cartoon illustrating the 
connection between tipf^ and tg: the electron exits the barrier with a negative velocity 
■v(ti,ph) = k(ti,ph) (directed towards the core). Its velocity gradually decreases and 
becomes zero at the classical ionization time {v{tB) = k(ts) = 0) - the time of birth tg. 



This section shows that if we relax the return condition and neglect the imagi- 
nary displacement between ti and Re(ti), we can obtain the same re-collision energy 
for real-valued canonical momenta and for real-valued return times. We shall call 
this approach the photoelectron model since it allows one to incorporate standard 
strong-field ionization concepts in a natural manner. The ionization amplitude would 
then correspond to creating an electron with a real- valued canonical momentum, and 
the imaginary part of the action integral would only be accumulated between t^ and 
Re{ti). 

In the classical model, one assumes that the electron trajectory is launched at the 
real 'time of birth' tg with zero instantaneous velocity. The electron instantaneous 
momentum at tg can be written as k(ts) = p -I- A{tB) ~ 0, where the canonical 
momentum p is a constant of motion (neglecting the core potential). The link between 
tB and p, p = —A{tB), links f^ via [p + A(ti)]^ = — 2/p to the complex-valued 
ionization time ti. In particular, for a linearly polarized laser field we have [A{tipfi)~- 
^{'ts)]^ = —'2Ip- Note that this t^p/i is in general difi^erent from the ionization time 
ti introduced in the previous section, since now the electron canonical momentum is 
forced to be real. The notation ij p/j stresses that this ionization time corresponds to 
photoelectrons, i.e. to electrons with real canonical momenta. Figure [L7| shows the 
mapping between the time of birth and the complex time tipf,. 

The photoelectron exits the tunnelling barrier at the real time, Re(ti p^), and since 
Re(ij p/i) turns out to be smaller than tg, the electron velocity at Re(ij p/j) is directed 
towards the core. It gradually decreases until becoming equal to zero at tg. The 
difference between Re(ti p^) and tg is small near the peak of the oscillating electric 
field, but increases as the field approaches zero. While the times tg are always spread 
within one quarter-cycle, as in the classical model, the times Re(ti p/j) are limited to 
a shorter fraction of the quarter-cycle, see Fig. |1.7| 

We now turn to the classical return time i^. In the original classical model, it is 
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0.2 0.4 0.6 0.8 1 

Real time of return, units of laser cycle 



Figure 1.8 Left panel: Physical picture of high harmonic generation In coordinate space 
and the meaning of different times. The electron enters the barrier at a complex time U 
and exits the barrier at the real time t'- = Ret;. Its velocity goes through zero at a (later) 
time ts- At the moment t^ the electron returns to the position It had at tg, and at the 
moment t It returns to the origin. Right panel: The blue curve shows the electron return 
energy at the moment t in the Lewenstein model, while the red curve shows the electron 
return energy in the classical three-step model. 



defined by the condition 



[p + A{t)] dr : 



fin 

JtB 



[-A{tB)+A{r)]dr^O. 



However, since the electron is already ofl^set from the origin at tg. 



Az = 



[Air) - A{tB)] dr, 



(1.87) 



(1. 






it does not return to the origin at t/j, see Fig. |1.8| 

The energy E{tji) = [A(tfi) — A{tB)]^ fi in the classical model is shown in the 
right panel of Fig. |1.8| with the cut-ofi^ at 3.17 Up + Ip. This cut-off is lower than 
in the quantum treatment, precisely because the electron has not yet returned to the 
core. The extra 0.32 Ip in the quantum cut-off law, 3.17 Up + 1.32 Ip, is due to the 
extra energy accumulated by the electron while covering the extra distance AzM 

Can we improve these results if we allow the photoelectrons to travel a bit longer 
and allow them to return to the core? Why do not we continue to monitor the electron 
trajectory at times t > tji and register their energy at the time of return to the origin 
tr^ph, ignoring whatever imaginary displacement they might have? There is just one 
problem with this plan: not all trajectories return to the core since we have limited 
the canonical momentum pp^ = —A{tB) to be no more than Aq. With this in mind, 
we shall take the energy at the closest approach to the origin as the return energy. We 
shall call this an improved three-step model or the photoelectron model. 



6) Interestingly, if one defines the experimental cut-off using the classical model, then the classical time- 
energy mapping is very similar to the quantum: Ir is very close to the real part of tj.. Since in the 
experiment the intensity is rarely known exactly, it is very difficult to differentiate between the classical 
(red) and the quantum ( blue) return energies in Fig. 1 1.81 
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Figure 1.9 Left panel: Energy of return for the Lewenstein model (red) and for the 
photoelectron model (green) vs real return time for /p = 15.6 eV, / = 1.3 ■ 10^* W/cm^, 
Yiuj = 1.5 eV. Right panel: Real part of the electron displacement in the photoelectron 
model. 



The model implies the neglect of the imaginary displacement and the minimisation 
of the real displacement between f , p/j and t^^ph- The imaginary displacement has 
to be neglected since we do not have imaginary canonical momenta and imaginary 
return times to cover for it. 

The photon energy resulting from the photoelectron model is Ef^^c (*) +^p = ^if > 
where S{pph, tr,ph^ ^i,ph) is given by Eq. (|1.21|l. It is in excellent agreement with the 



quantum photon energy (see Fig. |1.9[ left panel) for all those trajectories for which 
the real part of the electron displacement from the origin passes through zero. This 
is the case for the long trajectories and for most of the short trajectories, except for 
the shortest ones. These latter ones are 'born' at the end of the ionization window 
and contribute to the lowest harmonics, just above the ionization threshold. 

For short trajectories, the electron is decelerated by the laser field while returning 
to the core. Therefore, it needs a sufficiently high drift momentum to reach the origin. 
Since we have limited the canonical momentum p — -~A{tB) below A^, the shortest 
trajectories cannot quite make it to the core. For them, the time tr,p/i corresponds 
to the closest approach to the core. A non-zero real displacement yields a deviation 
of the approximate action S{ppfi,tr^pfi,ti^ph) from the real part of the exact action 
defined in the previous section, see the right panel of Fig. [L9] 

The action in this model is reproduced very well, since it is the time integral from 
the photon energy. Once the electron return energy is well-reproduced, so is the 
action, even if the end points t^, t[. are shifted. 

From the mathematical perspective, the photoelectron model implies that when we 
perform the integrals, we expand the action not at the exact saddle point, but in its 
vicinity. In particular, we shift the center of the expansion for the canonical momen- 
tum from the complex plane to the real axis. The error introduced in the integral by 
this procedure is minimized if the new expansion point lies within the saddle point 
region of the exact complex saddle point of the multi-dimensional integral. Thus, the 
difference Ap = pq — pp^ between the stationary point solution for quantum orbits 
Pq and the canonical momentum in the improved three-step model pp/j should be less 
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Figure 1.10 Left panel: Applicability region of tfie photoelectron model. The condition 
|A2:| < {tr — ti)^/^ specifies the region of return times (filled) where the photoelectron 
model can be used. The calculation is shown for Ip = 15.6 eV, / = 1.3 ■ 10^** W/cm^, 
huj = 1.5 eV. Right panel: Real and imaginary ionization times for the Lewenstein model 
(red), the photoelectron model with canonical momentum less than Aq (green), and for the 
canonical momentum not limited by this condition (blue). 
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than the size of the stationary point region: |Ap| < \d'^S/dp'^\^^''^ = {U — ti)^''^. 
We can estimate |Ap| as |Ap| = \Az/{tr — ti)\, where \Az\ includes the neglected 
imaginary displacement. This estimate yields \Az\ < (tr — ti)^''^. 

The left panel of Fig. |1.10| illustrates this condition for typical experimental pa- 
rameters (oj = 0.057 a.u., Ip = 15.6 eV, / = 1.3 • lO^'' W/cm^): the improved 
three step model cannot be applied for very short trajectories returning earlier than 
u)tj.^ph — 0.36 or for harmonics lower than N — 11. Thus, for this particular set of 
parameters, all above threshold harmonics are within the applicability conditions of 
the improved three step model. 

The right panel in Fig. |1.10| compares the ionization times resulting from the 
Lewenstein model and the photoelectron model of HHG. The ionization times coin- 
cide for the long trajectories. In this sense, the long trajectories indeed correspond to 
photoelectrons. The difference between the ionization times for the short trajectories 
is associated with the presence of imaginary canonical momenta in the Lewenstein 
model. For the shortest trajectories, the difference in the real ionization times is 
about 100 asec for the chosen laser parameters: the ionization window is wider for 
the photoelectron model. As for the imaginary component of the ionization times, 
they are smaller in the photoelectron model. Therefore, short trajectories are less 
suppressed in this model than in the full Lewenstein model. 

Mathematically, implementing the photoelectron model requires only one approx- 
imation - relaxing the return condition. Note that the requirement of perfect return to 
the origin is an artefact of neglecting the size of the ground state in the saddle point 
analysis. If we take into account the size of the ground state, then the return condition 
will naturally be relaxed: to be able to recombine, the electron has to return to the 
core within the size of the ground state. From this perspective, the extension of the 
Lewenstein model to real systems, including molecules, should go hand in hand with 
relaxing the return condition for its real part. 
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Measurement of ionization times might allow one to differentiate between these 
two models and to pin down the nature of the electron trajectories responsible for 
HHG. In particular, the interesting question is wether the complex momenta are the 
artefac of the 5-Uke initial state, or are indeed relevant for realistic systems. 

The accuracy of the first measurement [Shafir era/.| ( f2012) was sufficient to distin- 
guish between tg and i'j, but not high enough to distinguish delays between ^ (red 
curve. Fig. 1.10 left panel) and t[pf^ (green curve. Fig. 1.10 , left panel). 



1.8 

The multi-channel model of HHG: Tackling the multi-electron systems. 

In multielectron systems, there are multiple ways of energy sharing between the lib- 
erated electron and the ion. The ion can be left in its ground or in one of its excited 
electronic states. These options are referred to as different ionization channels. Mul- 
tiple ionization channels lead to multiple HHG channels: the returning electron can 
recombine with the ion in its ground or in one of its excited states. 

Multiple HHG channels present different pathways connecting the same initial and 
final state - the ground state of the neutral system - via different intermediate elec- 
tronic states of the ion. Thus, high harmonic emission in multielectron systems results 
from multichannel interference, see |Smirnova et al.\ ([2009a"), i.e. the interference 
of the harmonic light emitted in each channel. This interference naturally records 
multielectron dynamics excited upon ionization and probed by recombination, see 
[Smirnova et Q/.|p009a[ l. How important are these multiple channels? How hard is it 
to excite the ion during strong-field ionization? 

Strong-field ionization is exponentially sensitive to the ionization potential Ip, sug- 
gesting that after ionization the molecular ion is typically left in its ground electronic 
state. In the Hartree-Fock picture, this corresponds to electron removal from the high- 
est occupied molecular orbital (HOMO). However, multiple ionization channels can 
be very important in molecules due to the geometry of the molecular orbitals and the 
proximity of the excited electronic states in the ion to the ground state. 

The formalism described above, in the sections 1.1-1.7, is essentially a single- 
channel picture of HHG. It can be extended to multiple channels. 

First, we introduce the Hamiltonian of an N-electron neutral molecule interacting 
with a laser field: 



^c"-EE 



m i=l 



R-n 



N 



^e^=E 
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^i^ = -^F(i).d, = ^F(t).r,. (1.89) 

i i 

Here, the nuclei are frozen at their equilibrium positions Rm, the index m enumer- 
ates the nuclei with charges Qm, the superscript A'' indicates the number of electrons 
involved, T^ is the electron kinetic energy operator, V^ describes the Coulomb 
potential of the nuclei, Vje describes the electron-electron interaction, and Vf^ de- 
scribes the interaction with the laser field. Hats on top of operators are omitted. 

We will also use the Hamiltonian of the ion in the laser field, H^ ' , and the 
Hamiltonian of an electron interacting with the laser field, the nuclei, and the (TV — 1 ) 
electrons of the ion, H" = H^ - H^^-^\ 

The Schrodinger equation for the N-electron wavefunction of the molecule, initially 
in its ground electronic state *g (r), is 

*^(r,t = fo) = *^(r). (1.90) 

Similar to the single-electron case, its exact solution can be written as 

1*"^^) = -^ / di't/^(t,f')vf (t')<(i',to)*^(r) + <(t,io)l<>. (1.91) 

Jto 

Here the Uq and U are the N-electron propagators. The former is determined by 

idU^/dt = H^U^, (1.92) 

t/^(to,io) = l, (1.93) 

where H^ is the field-free Hamiltonian of the molecule: H^ = H — V^. The 
latter is the full propagator determined by idU^ /dt = H^U^ . 
The harmonic dipole reads 

Dit) = -z(f/^(t,to)<(r)!d| / df't/^(t,i')x 

Jto 

xVr{t')U^{t',to)^^{r)) + c.c. (1.94) 

Just as in the one-electron case (Eq. |1.13[ l, propagation without the laser field is 
simple as long as the energy Eg and the wavefunction of the initial state of the neutral 
molecule or atom are known: 

t/^(i',to)*f (r) = e-^^<'(*'-*°)*^(r). (1.95) 

Finding the full propagator U {t,t') is just as hard as solving the multi-electron 
TDSE. 

To simplify the analysis, we will make the following two approximations. First, 
we shall neglect the correlations between the electrons in the ion and the liberated 
electron after ionization. In this case, the full propagator factorizes into two in- 
dependent parts describing the evolution of the continuum electron and the evolu- 
tion of the ion in the laser field between ionization and recombination: U {t, t') ~ 
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f/*- ' (i, t')W(t, t'). Second, we will derive the results for short range potentials, 
just like we did in the single electron case considered above: U'^{t, t') ~ Uy{t, t'), 
and supply the corrections due to Coulomb effects. 

One can improve upon these two approximations by including the electron-electron 
correlations during ionization perturbatively, [Walters and Smirnova| ( |2010) , and by 
using the eikonal-Volkov states, [Smirnova et a/. | ([2008'), for the continuum electron, 
instead of the plane wave Volkov states. The eikonal-Volkov states include the laser 
field fully, the interaction of the continuum electron with the core in the eikonal ap- 
proximation, and also take into account the interplay between these two interactions 
(the so-called Coulomb-laser coupling, , Smirnova et al.\ ( (lQQl&) ). 

A consistent approach, which includes both electron-electron correlations and 
long-range effects in strong field ionization can be developed within the time- 
dependent analytical R-matrix (ARM) method ( [Torlina e?a/.| ( f2012j ). This method (i) 
splits the configuration space into the inner and outer region, uses quantum chemistry 
in the inner region, (ii) the eikonal-Volkov propagation in the outer region,( [Smirnova| 
\et a/.| ( |2008j ), and (iii) the Bloch operator CBloch (19571) to match the solutions in 
two regions. 

Moreover, if we can factorize the dipole response into the usual steps - ionization, 
propagation, recombination, we can think of improving each of the three steps sepa- 
rately, e.g. by using improved ionization and recombination amplitudes that include 
the electron-electron correlation beyond the perturbation theory. 

Just like in the one-electron formalism considered above, we will introduce the 
identity resolved on the momentum states of the continuum electron, but now we 
also have to include the electronic states of the ionP^I 

1= /"dp^Alrj(^-i)®pr){n(^-i'®pr|A, (1.96) 

where A denotes the antisymmetrizing operator. 
The harmonic dipole becomes 

D(i) = -i{*fldl /"*dt'e^^«(*-*') /"dpt/(^-i)(t,t')|n('^"'^)x 

X t/f.(t,t')|p?)(prm(^-l)il//'(t')l<) + C.C. (1.97) 

Note a crucial change compared to the single-channel case (Eq.fTTSj: the appear- 
ance of the laser-induced dynamics between the bound states of the ion, described by 
the propagator U^"^^> (i, i')|n^ '). This dynamics can be calculated if the dipole 
couplings, dmn, between all essential states, as well as their eigenenergies. En, are 
known. 

Consider, for example, the case of an N2 mo lecule with three essential states in the 
NJ ion, denoted as X, A and B, see Fig 1.11 The time-dependent transition ampli- 



tudes amnit, t') between the state rS^ ^populated at the moment t' and the state 

7) Here we use the field- free states of the ion. If the limited amount of basis states is used, then one should 
try to find the optimal "laser-dressed" basis. 



Thomas Schultz and Marc Vrakking: Attosecond and Free electron Laser Science 
Chap. 1 — 2013/4/10 — 2:48 — page 33 



Populations of ionic states 
|X|2 




Am, 



1.3 "eVf 



X2X+ 



3.5eV 



Ip~15.6eV 



Time, laser cycle 

Figure 1.11 Left panel: Sub-cycle dynamics in the nJ Ion aligned at 9 = 50° to the laser 
field polarization: populations of the field-free Ionic states X (blue), A (red), and B (green) 
in a / = 0.8 • lO^"* W/cm^, 800 nm laser field. Right panel: Electronic states of the N J Ion. 
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(N-1) 



at the moment i are given by ar„n(i,i') = {m(^-i)|t/(^-i)(t, t')!"-*^^"^)). 



It is a solution of the following system of differential equations: 

rf(A("')_r„, w.M.(n) 



dt 



[H + V(i)]A 



(1.98) 



^Ei 



\ 
where, for our three ionic states, the Hamiltonian of the ion is H = E2 , 

V E3) 
with the energies En of the three states. The interaction between these three states is 
described by the matrix of the laser-induced couplings, Vmn{t) = — dmn ■ F(i), that 

/ Vi2it) Vi3{t)\ fain{t,t')\ 

is V(f) = V2i(t) V23(i) . Finally, A^") = a2n (*,*') is the vector 

\V31it) V32(i) / V3„(i,t')/ 

describing the population amplitudes of all essential ionic states, starting from the 
state n(^"^) at time i'. 

Let us introduce channel specific Dyson orbitals ^^^(r) = {n'^^^^|*^(r)). 
These are the overlaps between the A'^-electron wavefunction of the ground state of 
the neutral and the (TV— l)-electron wavefunction of the ionic state \n^ '). Let us 



assume that the dipole operator that starts ionization at the moment t' in Eq. (1.97 1 
acts only on the electron, that will be liberated (i.e. we neglect the exchange-like 
effects in ionization). In this case, the multielectron dipole Dmn, which corresponds 
to leaving the ion in the state n^ ' after ionization and then recombination with the 
ion in the state m^ ' , can be re-written in a form very similar to the one-electron 
case (Eq. ^L24\ y. 



D 



(mn) 



(i) 



i f dt' /°dpd;,(p + A(t))a™„(t,t')e^'^"'P'*'*'^ 



xF(i')d„(p + A(i')), 



(1.99) 
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d„(p + A(i)) = {p + A(f)|dl*,?), 

d„(p + A(i)) = (p + A(t)|(n(^-l'ld|*f >, 

1 /■* 

Sn{p,t,t') = - [^ + A(r)]2dr + /p,„(i - t'). 

This expression is remarkably similar to one-electron dipole (fL24l. The transforma- 
tion similar to (|1.22^ is also valid in this case, yielding 



D('"")(i) = i 



j dt' fdpdUp + Mi))amn{t,t')e~'^^P^'^'"> X 



xT„(p + A(t')), (1.100) 



T„(p + A(i)) = 



l2 



7, r lp,n 



(p + A(i)l*^), 



where Ip^n is the ionization potential to the state n of the ion and the matrix amn (t, i') 
is calculated while setting En to zero. 

The total harmonic signal results from the coherent superposition of the dipoles 
Dmn associated with each ionization-recombination channel: 

D(i) = ^D('"")(i). (1.101) 

m,n 

Substantial sub-cycle transitions, such as those shown in Fig. |l.ll| for the N^ ion 
in typical experimental conditions, have a crucial impact on the harmonic radiation. 
They lead to the appearan ce of th e cross-channels in HHG (the off-diagonal elements 



j-)(mn) fojmi^n in Eq ( 1.101 ) since the state of the ion changes between the ion- 
ization and the recombination, see Fig. |1.12[ These channels are indeed substantial 
in high harmonic generation from the N2 molecules, see |Mairesse et i2/.] ( |2010[ l, as 
illustrated in Fig. |l.ll| 

In the recent literature on high harmonic generation one can often come across a 
rather loose language, which refers to different ionization and recombination channels 
as associated with different Hartree-Fock molecular orbitals. This language should 
not be taken literally as a statement on the applicability of the Hartree-Fock picture 
and on the physical reality of the Hartree-Fock orbitals as observable physical quan- 
tities. Loosely speaking, removing an electron from the highest occupied molecular 
orbital (HOMO) creates the ion in the ground state. Removing an electron from one 
of the lower lying orbitals (e.g. HOMO-1, HOMO-2) creates the ion in one of its 
excited states. Thus, the reference to the orbitals should only be understood as a lan- 
guage for describing ionization and recombination channels associated with different 
multielectron states of the ion - and those are physically relevant and observable. In 
the orbital language, electron removal from an orbital creates a hole in this orbital. 
The laser induced dynamics in the ion, moves the hole between the orbitals in the 
time window between ionization and recombination, see Fig. |1.12| 

Application of the saddle point method in each channel leads to the following half- 
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Ion Channel XA 
\_/ ionization 





recombination 



Ion* Channel KA 
-I i— ionization 




recombination 



Figure 1.12 Left panel: Cross-channel in HHG associated with ionization from and 
recombination to different orbitals. This channel is due to real excitations induced by the 
laser field between ionization and recombination. Right panel: Diagonal channel in HHG 
associated with ionization from and recombination to the same orbital. 



cycle dipole for the given ionization - recombination channel: 



OionVPs, ti) — 



27r 



nl/2 



iSips.t/'^t'r") 



T„(ps 



^prop 



{Ps,t,ti) = 



(27r 



,3/2 



iS(p,,t,<«')^ 



(i-i»)3/2 
d*r„{ps+A{t)). 



i{t,t. 



A(tp^)), 



(j)x 



(1.102) 
(1.103) 

(1.104) 
(1.105) 



Neglecting the sub-cycle dynamics in the prefactor of a"Q„(ps, f^) we can substitute 
Eq. ( |1.103[ l by the following expression, which includes the Coulomb effects in ion- 
ization (the sub-cycle Coulomb effects see|Torlina and Smirnova|(|2012)): 



i{Ps,k) = 2 



2k' 
~F 



-'yuj 



^~iS{ps.t'i,U) 



-'-\/l+~- 



T„(ps+A(ii)).(1.106) 



A term similar to Tn also arises within the time-dependent analytical R-matrix ap- 
proach applied to multi-channel strong field ionization |Torlina et al. | ( [2012^ . However 
in |Torlina et a/.|pbl2[ l, the radial integration is removed due to the use of the Bloch 
operator. Thus in |Torlina et fl/.|p012) , the pole in T„ does not arise even when the 
long-range potential is taken into account. Function Tn (ps + A(f ^ ) ) encodes the an- 
gular structure of the Dyson orbital in the asymptotic region, which is more complex 
than the one arising in the asymptotic of the atomic wave-function ( |1.56[ l leading to 
Eq. ( |1.58^ . The simple expressions for the asymptotic of the Dyson orbital for small 
molecules can be found in Murray et a/.|p01 1) . 

Here, we have considered the harmonic dipole on the real time axis. Note that 
the propagation amplitude is modified compared to the one in the one-electron case 



Thomas Schultz and Marc Vrakking: Attosecond and Free electron Laser Science 
Chap. 1 — 2013/4/10 — 2:48 — page 36 

36 I 



(Eq. 1.72[ to include the laser - induced dynamics in the ion amn (t) ■ The full dipole 



for each ionization-recombination channel is the sum over the different half-cycles 
and the harmonic spectrum results from the FFT of the full dipole DmniNco): 

D('"")(i) = ^D(j'"")(t), (1.107) 

j 



: / dfD 



D^'""^(7Va;)= / dfD^'""^(f)e"'"^ (1.108) 



The complete harmonic response is obtained by adding coherently the contributions 
of all ionization - recombination channels. 



1.9 
Outlook 

Having factorized the dipole, we can use improved amplitudes for each step. These 
are the key components of the current theoretical work in high harmonic spectroscopy 
of molecules. 

Improving ionization. Improved ionization amplitudes can be taken from semi- 
analytical and/or numerical approaches. The task is to define the function T(p + 



A{t)) for a realistic system and include long-range and polarization ( Sukiasyan et al. 



l |2010 )) efl^ects. For example, one can use the results of |Murraye? al. \ ( |20 11) , where 
the ionization amplitude is represented as: 

«ro„(Ps,t,) =R,„(/p,i^)e-*^(P-*'^*'). (1.109) 

The exponent describes the sub-cycle dynamics of strong-field ionization, i.e. is the 
same as for the atomic case and a short-range potential. The prefactor R|m(^p, F) 
accounts for the influence of the core potential and the shape of the initial state on the 
ionization rate. For atoms, this prefactor has been derived in the seminal papers of 
Perelomov, Popov and Terent'ev (see [Perelomov et al] \ 1 966{ 1967^; Perelomov and| 



|Popov| (1967) ; [Popov et a/.| ( |1968) ) and improved in [Popruzhenko et a/.] ( |2008^ . A 



simple recipe for incorporating their results into the sub-cycle ionization amplitudes 
can be found in |Yudin and Ivanov| ( |2001) . Fully consistent treatment of long-range, 
polarization and orbital effects can be developed within the time-dependent analytical 
R-matrix (ARM) approach, see |Torlina and Smirnova| ( |2012) ; [Tbrlina e?a/T| ( |2012^ . 

Improving propagation. In addition to the dynamics in the molecular ion, in- 
cluding the laser-induced transitions between different ionic states, the second most 
important modification of the propagation amplitudes is the incorporation of possi- 
ble transverse nodal structure in the continuum wavepackets. The nodal planes in 
the continuum wavepacket arise during tunnelling from bound states. For example, 
in the CO2 molecule, the HOMO and the corresponding Dyson orbital have nodal 
planes parallel and perpendicular to the molecular axis. Consequently, in the case of 
tunnel ionization with the molecular axis aligned parallel to the polarization of the 
ionizing field, the nodal plane will not only reduce the ionization rate, but will also 
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be imprinted on the shape of the electronic wavepacket that emerges after ionization. 
Propagation between ionization and recombination will lead to the spreading of the 
wavepacket, but it will not remove the presence of the node as the wavepacket re- 
turns to the core |Smirnova et Qr| ( |2009a,c,bj . Clearly, this aspect of propagation is 
important for the recombination amplitude. 

Consider, for example, ionization from a state with angular momentum L — 1. 
Its projection on the laser polarization is either Lz — Q (no nodal plane along the 
electric field) or Lz = 1 (nodal plane along the electric field). After tunnelling, in 
the plane orthogonal to the laser polarization, in t he momentum space the continuum 
wavepackets are proportional to (^) , (see Eq 



1.62 



1.66i: 



*L,=o(p_l) oc e" 

^L^=l{Pl_)^—e^'e-'^\ (1.110) 



K 



where r — Im(ii), k = \/2Ip, and 4>p is the angle between p_\_ and the a;-axis. 
As we can see, the presence of the nodal plane for Lz — 1 leads to the additional 
term pj^/k. We now propagate these wavepackets until the recombination time tr. 
Fourier transforming back into the coordinate space, in the plane orthogonal to the 
laser polarization, we get 

2 

^L.=l(p)« P e"^e~^(^ (1.111) 

where p is the transverse radial coordinate and 4> is the angle between the radial vector 
and the x-axis. Recalling that x = pcos cj>, we see that if we combine the Lz = ±1 
states to form the real- valued spherical harmonic px, the presence of the nodal plane 



efi^ectively changes the dipole operator d to d ■ x/{K{tr — ti)). In Smirnova et al. 
( |2009a|c|b) such modifications of recombination operators has been used to account 
for the appearance of nodal planes. 

In most experiments with molecular HHG to-date, the alignment distribution is 
rather broad. Even if the molecular ensemble is, on average, aligned parallel to the 
laser polarization, for most molecules the characteristic alignment angle would be 
sufficiently different from that associated with the nodal plane. In this case, the rela- 
tive importance of the nodal planes in recombination is reduced. However, for well- 
aligned molecular ensembles this would become a significant factor. 

Improving recombination. The recombination step can be significantly improved 
beyond the SEA, if one uses the recombination dipoles dm(ps + A(t)) calculated 
using ab-initio approaches. Eor example, the quantitative rescattering theory (see 
|Lin et aL\^20lO) ) relies on using the Schwinger variational method to calculate the 
field-free recombination matrix elements. Alternatively, one can use the R-matrix ap- 
proach (see parvey and J. |P009[ | ; [Harvey et aL\\20l2\ ) . Both allow one to incorporate 
the full complexity of the recombination process, including the channel coupling due 
to the electron-electron correlation and automatically include the exchange effects in 
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recombination [Gordon and Santra| ( |2006) ; [Patchkovskii ef a/.| ( [2006| l; ?. The draw- 
back of these methods, at the moment, is the absence of the laser field in the calcula- 
tions of the recombination amplitudes. This approximation breaks down in case of a 
structured continuum |Patchkovskii et gT] ( |2012j , common for many molecules. The 
impact of the IR field on such continuum states has been recently demonstrated exper- 
imentally 'DtteFari ([20131, substantiating the prediction of |Patchkovskii et ar\\2Q\2) . 
In the approach described by |Smirnova et al.^ l ^2009aj , the eikonal-Volkov approxi- 
mation for the continuum states was used to obtain improved dipoles in the single- 
channel approximation with exchange. The eikonal-Volkov approximation fully in- 
cludes the interaction of the continuum electron with the laser field, but the inter- 
action with the core potential is only included in the eikonal approximation, and 
the correlation-induced channel coupling is neglected. Improving the recombination 
amplitudes to account for all these effects - the channel coupling due to the electron- 
electron correlation, the core potential, and the laser field, is one of the key theoretical 
challenges today. 

With each of the three steps in the harmonic response improved, the original SFA- 
based theory turns from purely qualitative into a little more realistic. The separation 
of the three steps, crucial for our ability to improve each of them separately, benefits 
from the high intensity of the driving field and the large oscillation amplitude of the 
active electron. The high field intensity also lies at the heart of the main difficulties 
in building an adequate theoretical description. Nevertheless, the effort is worth the 
investment: the combination of attosecond temporal and Angstrom spatial resolution 
is extremely valuable. High harmonic spectroscopy appears to be well suited for 
tracking the multielectron dynamics induced by the ionization process. 

It is very attractive to replace the ionization step induced by the IR field with the 
one-photon ionization induced by a controlled attosecond XUV pulse, phase-locked 
to the strong IR field (see |Schafer et al. |p004[ l). The latter would drive the continuum 
electron. Such an arrangement should allow one to move from dealing with outer va- 
lence electrons to dealing with inner valence and deeper lying electrons. This appears 
to be an exciting regime for tracking the hole dynamics (Lun nemann ef a/.| ( |200 8^) ini- 
tiated by inner- valence or deeper ionization. Importantly, for deeply bound orbitals, 
the effect of the IR driving field on the core-rearrangement and the hole dynamics 
should be substantially less than for the outer-valence electrons. 

High harmonic spectroscopy helps to record the relative phases between different 
ionization channels, which provide information about electron rearrangement during 
ionization and define the initial conditions for the hole migration both in the inner 
valence l |Lunnemann et al. | ( |2008 1) and outer valence ( |Smirnova et al. |p009a[ l) shells. 
These initial conditions are sensitive to the frequency, intensity and duration of the 
ionizing pulse, opening opportunities for controlling hole migration and, possibly, 
related chemical reactions [Weinkauf et al.\^991\ . 
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1.11 



Appendix A: Supplementary derivations 

In this Section we prove that the transformation ( |1.22[ l 

r-t 

I dt'e"*^(P'*'*)F(f')d(P + A(i')) = 

Jto 

ft 

= I di'e-'^(P'*'*)T(p + A(f')), 

Jto 

(P + A(i'))' 



T(p + A(i')) 



+ h 



(p + A(t')|5>, 



is applicable in case of high harmonic generation. By definition 

rt 



I di'e-'^(P'*'*')F(f')d(p + A(i')) 

Jto 



(1.112) 



(1.113) 



(1.114) 



''-' / dt'{^l(t'-t)\-VL\g[t')). 



Adding and subtracting the kinetic energy operator ^ (Becker et al. 
|akin and Kuchiev| ( |1997[ l) we obtain: 

= / 

Jt, 



1 2002b I, 



Grib- 



^* S^ 7m \9{t)) 

'to 



'^^2^ 



P + A(i')) 



.(S*r(t';t)l 
Here we have used that —i yp - 

the first term in Eq.l|1.116| we obtain: 



;*p(t';i)l5(i')) 



(1.115) 
(1.116) 



T + ^lI {*p (f ; t)\- Integrating by parts 



/' 

Jto 



dti ^, \g{t)) 



dt' 



(1.117) 
(1.118) 
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The transformation (1.1 12 i can be recovered using i ^y^, " = —Ip\g{t')) and taking 



"3F 



/ 



into account that the boundary term (*p (t ;t)\g(t ))\l^-^ does not contribute to the 
high harmonic dipole. Indeed the contribution of the boundary term to high harmonic 
dipole is: 

dpd*(p + A(i)) [(*^(t;t)l5(t))-(*^(to;t)l5(to))l = (1.119) 

dp{git)\d\^^{t;t)){^^{t;t)\g{t)) - (1.120) 

dp{g{t)\d\^^{t-t)){^^{to;tMto)}. (1.121) 

The term ( |1.120| ) is equal to zero, the term ( |1.121| ) tends to zero when to — > — oo. 
Indeed, 

dp(.g(t)ldl*^(t;t))(*^(t;f)l3(f)) = {g{t)\d\g{t)) = 0, (1.122) 



while the second term l |1.121| ) is: 



dp(g(i)id!*^(t;t)){*^(0;i)l5(0)) = (1.123) 



/ 



dpe 



-*i /4 dr[p+A(r)]^^^^^^|^|p ^ A{t)){p\g). (1.124) 



This term corresponds to the projection of the ground state onto the basis of plane 
waves at to — >■ — oo followed by recombination of the resulting oscillating wave- 
packet back to the ground state at time t. Spreading of the free electron wave-packet 
over infinite time {t — to) — > +cxd makes this projection negligible. 



1.12 

Appendix B: The saddle point metliod 

The saddle point method is one of the key techniques in the analytical strong-field 
theory. It is an asymptotic method, which allows one to analytically evaluate the 
integrals from highly oscillating functions, such as the integral in Eq. (fTTT9}. 

1.12.1 

Integrals on the real axis 

How would one calculate the following integral, 

(•b 

f{x)e^^^'='>dx (1.125) 

for some smooth functions f{x) and h{x), without knowing much about them, or if 
they look ugly and complicated? All we know is that they are real-valued functions 
on the real axis x. 
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In general, one could think that there is not much one can do. Fortunately, this 
is not the case if the positive and real A is large, A 2> 1 - then, the integral can be 
calculated. 
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1.12.1.1 Contribution of the end points 

The first idea that comes to my mind when looking at such an integral is to try in- 
tegration by parts. This approach works just fine under certain circumstances, see 
below. The first stumbling block meets you right at the gate: how does one integrate 
by parts if both h{x) and f{x) are unknown? 
The trick is simple: 



I = 



/' 

J a 



fi^)e 



Xh{x) 



dx 



'"^Jgy^'-^"""" 



1 f{x) xh(x),b _ 1 /■ ,XHx) \± ( f{x) 
Xh'(x) ''' A./^ [dx \h'{x) 



(1.126) 



We started with an integral that did not have a small parameter 1/A in front. Now we 
have two terms: the first comes from the contributions at the end points. The second 
term is another integral, now with a small parameter in front. Dealing with it in the 
same way as with the original integral, we will get terms proportional to 1/A^, and 
so on. 

Thus, we conclude that the main contribution to the integral comes from the end 
points, and is given by the first term: 



I = 



f{x)e^''''-''Ux 

f{b) \h(b) _ fja) \h{a) 



h'(h) 



h'{a) 



+ 0(A-2). 



(1.127) 



This result is applicable unless there is a problem with the second term in Eq. ( |1.126p 
- the integral 



dxe 



\h{x) 



d f fix) 



dx \h'{x) 



(1.128) 



The problem arises if h'{x) = somewhere between the two end points of the inte- 
gral. What do we do then? Obviously, the points where [f{x)/h'{x)] diverges can 
bring major contributions to the integral. 

Given that A » 1, the way the exponential function changes between a and b is 
most important. The first possibility is /i' / in the integration interval. Then, 
the integral is accumulated at the end points, and the end point where h{x) is larger 
dominates. In general, for an exponential function e '^^ the main contribution to 
the integral will come from the region where it reaches its maximum value - and 
hence where h'{x) — 0. 
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Figure 1.13 Integral from a bell-shaped curve with a maximum at xq and a width of Ax. 

Suppose that, somewhere between a and 6, the derivative h' — 0. If the function 
h{x) has a minimum, the contribution of this minimum won't be competitive with 
the contributions from the end points (remember that A is large and positive). But if 
it has a maximum, then the main contribution to the integral comes from the region 
near the maximum. The way to handle this situation is described in the next section. 



1.12.1.2 The Laplace Method 

Let us consider an integral from a function f{x) shown in Fig. 1.13 The function is 



bell-shaped, has a maximum at the point xq, where its first derivative is, of course, 
equal to zero, and quickly falls off to each side of xq. 

Calculation of this integral is very simple - all we need is to find the effective width 
Ax of the bell-shaped curve, and then the integral is 



fix)dx = f{xo)Ax. 



(1.129) 



Let us first try some simple estimates of the width Aa;. In order to do it, we expand 
f{x) around xq in a Taylor series, remembering that the first derivative is zero at this 
point: 

fix) ^ f{xo) + lf"{xo)ix - xof = f{xo) - l\fixo)\{x - xof. (1.130) 

Notice that I have explicitly used the fact that the second derivative at the local max- 
imum is negative. 

A potential candidate for the width Ax is the full width at half maximum (FWHM). 
The half-width Ax/2 at each side is given by 



fixo) -l\.f{xo)\{Ax/2f ^ fixo)/2. 



(1.131) 
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This gives us Ax = a/^/o/I/"!- 

A more accurate calculation of the required width comes from the following trick, 
which will also smoothly bring us into the saddle point method 

7= [ f{x)dx= /e'"-^(=")da;. (1.132) 



This transformation allows us to reduce the integral to a familiar Gaussian form. We 
proceed by expanding ln/(a;) in a Taylor series, remembering that /'(xq) — and 
fixo) = -\r{xo)\: 

.i"/(-)d.^ /e'-^-^^^^-^^'^^dx. (1.133) 



Recalling that the integral from a Gaussian is 

e'^'^'da; = a/^, (1.134) 

and setting the limits of our integral to ±00, we get the final answer 



/- 



00 



I = f{xo)^27vf{xo)/\r{xo)\. (1.135) 

As you can see, the width Aa:: turned out to be pretty close to the FWHM. 

1 .12.1 .3 Saddle point method: the steepest descent in a complex plane 

We now move to the saddle point method which is used for integrals of complex- 
valued functions: 

e^^^'^Uz. (1.136) 

where A is large and positive, and the rest is hidden in f{z). The integral is to be 
taken over a contour C, and the only good thing about this contour is that its ends, 
somewhere far away from the place of action, do not contribute to the value of the 
integral. 

There are assumed to be no poles, so that we are allowed to deform the contour C 
as we wish. The key of the steepest descent is a clever modification of the integration 
contour. 

First, note that a complex function f{z) has a real part and an imaginary part, 
/(z) = u(z)+iv{z) = u{x, y)+iv{x, y), where x and y are the real and the imaginary 
parts of z, z = X -\- iy. 

Let us now look at the integral more closely and recall the previous section, where 
the integration was based on expanding the function around a maximum and reducing 
the integrand to a Gaussian. In our case we have a function oxp( Au+iAw) that changes 
its absolute value very rapidly due to the \u part. It also oscillates rapidly due to the 
\v part. The trick of the steepest descent is to modify the contour of integration in 
such a way that it will go through a place where the real part u climbs to a maximum 
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along the contour and then quickly falls, while the imaginary part v stays constant 
along the same contour, freezing any fast oscillations. 

It may not be obvious at first glance that such a modification of the contour is 
possible. But it is. 

We start in a manner entirely analogous to the previous section. Let us assume that 
the function f{z) has a zero derivative at the point zq — XQ,yQ, where xq and j/q ^e 
coordinates in the complex plane; the point zq lies somewhere between the left and 
the right ends of the contour C 

If fz{zo) = 0, then both the real and the imaginary parts of / must have zero 
derivatives there: 

Ux = Vx =0, Vy =Uy =0. (1.137) 

Thus, not only at zq the absolute magnitude of our function goes through an ex- 
tremum, but also the oscillating part is stationary. Another important observation is 
that the gradients of the two functions, Vu and Vn, are always orthogonal to each 
other: 

\/v ■ VU = UxVx + UyVy — 0. (1.138) 

This is the consequence of the Cauchy-Riemann conditions: 

Ux—Vy, Vx — —Uy. (1.139) 

The gradient points into the direction for which the function changes. If we move 
along the gradient of u, following the path of its steepest rise and fall through the point 
zq, we are also moving orthogonal to the gradient of v. Thus, v will stay constant, 
and fast oscillations are frozen. We see that the desired modification of the contour 
is indeed possible. 

How should the landscape of u{x, y) look like? Due to the same Cauchy-Riemann 
conditions the functions u and v can only have saddle points at zq: 

Uxx+Uyy=0, Vxx+Vyy=Q- (1.140) 

Real mountain peaks, which go down in all directions, only happen at singularities, 
and we decided that there would be no singularities in f{z). 

Therefore, the landscape of the function u around the point zq must look something 
like shown in Fig. |1.14| 

All we have to do now is to find the correct path of the steepest descent through the 
saddle point, such that u will rise as quickly as possible as we approach the saddle 
point and then decrease as quickly as possible as we leave the saddle point. The 
Cauchy-Riemann conditions promise us that, while we are at it, v will stay constant. 

Let us expand f{z) in a Taylor series around zq, remembering that /'{zq) = 0: 

f{z)^f{zQ) + ^f{zo){z^zof. (1.141) 
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y=Im z 




Landscape of 
u(x,y)=Re f(x) 



direction of 
descent 



x=Re z 



Figure 1.14 Saddle point method. The landscape of u{x,y) = Kcf(z) around the point zq 
where /' = 0. 



Of course, f"{zo) is a complex number, which we will denote as /"(zq) = a oxp{id). 
If our path traverses the saddle point at some angle (f>, then z — zq = p exp(i0) and 



1 r//f ^r \2 1 2 i(0 + 2<h) 

^f {zo){z-zo) = -ap e^ ^ f'. 

Now, the simple trick is to choose the angle (p properly - we set 

i(e+20) _ 1 



(1.142) 



(1.143) 



and keep the angle <j> given by the above condition fixed, changing only p, so that 
dz = d{z — Zq) = exp{i<j})dp. 
If we do this, the integral along such a path will look as 



7 = e 



^f(za)^i'P 



e 2 (ip^ 



(1.144) 



where the deformed contour C' is going through the saddle point as a straight line at 
an angle 4>. Note that, indeed, there are no oscillations along such path, and the real 
integrand decays as a Gaussian. 

The integration limits are now extended to plus and minus infinity and the integral 
is done: 



V Aq 



(1.145) 
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Recall that as \f"iza)\. 

At this point we are almost done, but three important remarks are still in order. 

First, there is ambiguity in the definition of the direction (f> from exp{i{d + 2(/>)) — 
— 1. Indeed, the total angle 6 + 2(j) could be both plus and minus tv. Thus, formally, 
without looking at the landscape shown in Fig. |1.14| we have a choice of two 0: 

(j>i = -e/2 + TT/2, (j>2 = -d/2-Tr/2. (1.146) 

The whole idea of the method is to choose such a direction that you never have to 
cross the 'mountains' in the landscapes of u{x, y) and v{x, y). You should choose the 
direction (deforming the contour C) that takes you from the valley, through the saddle 
point, into another valley. Otherwise, you will also have to include the contributions 
of the 'mountains' into the integral. Usually, it is the first choice that works, but 
one should take a look at the landscape and check. The wrong option will go in an 
obviously wrong way, crossing into the tops of the mountains rather than staying all 
the way in the valley and smoothly climbing to the saddle. We shall see an example 
of it in the next section. 

Second, if there are several saddle points, i.e. f{z) has many points where its 
derivative is zero, the integral will be the sum of the contributions from all these 
points. Then the individual phases (f) for each saddle point become very important. 

Third, there is a modification of the prefactor when dealing with multi-dimensional 
integrals: 

(1.147) 

(1.148) 

where fzz is the Hessian matrix (the matrix of the second derivatives of the function 
f). 

1.12.2 

Stationary phase method 

The stationary phase method is a simple application of the saddle point method to a 
function with a purely imaginary phase: 

/= j g{x)e'-^f^'''^dx, (1.149) 

where g{x) is a benign, very slowly changing function which does not do much - 
just makes sure that the integrand goes to zero at infinity. The constant A is again 
real and positive, the integral is supposed to be performed along the real axis and the 
function f(x) is purely real on the real axis. Intuitively, it is clear that if the function 
exp(JA/(a;)) is oscillating very quickly, its integral averages to zero unless there are 
some points where the oscillations freeze. These areas are the regions where the 



1 = 


Jc 

,p-\"/%A/(.„) 1 




^ ^ ^ V-fM^o) 
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stationary 
phase point 




Figure 1.15 Stationary phase method. The integral comes from the area where the 
integrand does not oscillate as much. 



phase of the oscillation, f{x), stays nearly constant, i.e. areas around the point where 
the derivative turns to zero, /' 



1.15 



0, see Fig. 

The problem can be turned into that studied in the previous section. Again, suppose 
that the derivative /' = at some point xq. We use the same Taylor expansion around 
this point and denote x ~ xq as, say, ^. The integral is approximated as 



I = g{xo)e 



iA/(. 



"7' 



i>^f"{^ 



-de 



(1.150) 



and we will assume that /" = q > (/ is a real- valued function and xq is on the real 
axis, hence /" is real). The case /" < is handled in an identical manner. 
Calculation of the integral 



/ 



e ^ at. 



(1.151) 



follows the exact prescription of the saddle point method. Obviously, the phase 6 of 
the second derivative (see previous section) is ^ = 7r/2 (i.e. ia — aexp(J7r/2)) and 
the contour of integration has to be turned at an angle 4> to the real axis, such that 
+ 2(1) = ±7r. This yields the two possible choices of 0: (j) — +-k/4, and <^ = — 37r/4; 
the answer for the integral is: 



g{xo)e 



(1.152) 



To find the correct choice of <^, one has to look at the landscape of u{x, y) = 
Re(iz^) — —2x1/. The landscape is shown in Fig. 1.16 The correct choice is ob- 
viously the first one, <^ = +7r/4, the second would mean that the contour has to be 
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Landscape of 
u(x,y)=Re iz-= -2xy 

(|)=7l/4 - correct 
direction of descent 



x=Re z 

-371/4 - wrong 
direction of descent 



Figure 1.16 Stationary phase method. Correct and incorrect paths of the steepest descent 

for f{z) = iz'^. 



deformed as shown in Fig. |1.16| with the dashed line, going through high mountains 
on the way to the saddle to cross it in the opposite direction oi 4> = — 37r/4. 
So, the final answer is 



I = 9{xo)e \i—e 



i\f (xo) J ^t^/4 



(1.153) 



1.13 

Appendix C: Treating tlie cut-off region: regularization of the divergent 

stationary phase solutions 



In this subsection we briefly outline the idea of the so-called uniform approximation - 
one of the approaches for handling the merging stationary points. The regularization 
involves two steps. First, we need to find a specific real return time tr = tro and the 
associated ti = tio, ps = Pso , such that d'^S{tro,tiQ,pso)/dt'^ = 0. In practice, one 
can simply pick the real return time corresponding to the cut-off energy. The next 
step requires the expansion of the total action in Eq. ( |1.2ip around t = i^o in a Taylor 
series up to the third order: 

S(t,t,o,Pso) = S{tro) + {t- tro)S'u + ^^^^^ S't't + i^^^s't'It, (1.154) 
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where all the derivatives of S{t, iiOiPso) are taken at trQ. Finally, one substitutes the 
expansion (fTTT54l into the harmonic dipole. 



D{Nuj) oc 
and uses the Airy function 



fllgi'^'^t^--iS{t,tio,Pso)+iNu]t _|_ g g 



r 



dt cos{at ± xt) 



(3a)i/3 



Ai 



±x 



(3a)i/3 



(1.155) 



(1.156) 



Now we introduce the 'cut-off harmonic number' A'^n and the distance from the cut-off 

N ^ N - Nq: 



NqLU = Erec{tro) + Ip, 



(1.157) 



here Erec(tro) = PsO + A{tro) is the re-collision energy at time tro, and A^o does 
not have to be integer. The dipole near the cut-off is expressed via the Airy function: 



flllJ-'^'^'t g-'>'S{t,tiO,Pao)+iNu:t : ^ ^ _ 



SO that 






ANuj 
i(x/2)i/3j 



d^cos{^f±ANuj^), (1.158) 



(1.159) 



where x = ~'S'"/t(iro) and can be estimated by x — v(tro)Fouj, given that 
Ftix) — ^01^ and F{x) = 0. Using the asymptotic expansion of the Airy func- 
tions, we obtain simple expressions for the dipole just before and after the cut-off. 
Before the cut-off of the harmonic spectra (for AA'^ < 0), the dipole oscillates, 
Ai ~ cos{-{ANuj)^^'^{8/9x)^^^), after the cut-off, the harmonic dipole expo- 
nentially decreases, Ai ~ exp(— (AA''a;)^"(8/9x)^ )■ The oscillations of the 
harmonic dipole before the cut-off are due to the interference of the short and the 
long trajectories. 
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1.14 

Appendix D: Finding saddle points for the Lewenstein model 

In section 1.5 we have described how one can find all saddle point solutions in 
the Lewenstein model for a fixed harmonic number A'^. Here we present an alter- 
native and equivalent approach of finding the saddle point solutions, i.e. solving 
the Eqs. ( |1.36|1.37|1.38^ , which can be used in all cases, but is particularly con- 
venient if the Fourier transform is performed numerically. The idea is to solve the 
Eqs. ( |1.36|1.37|1.38^ 'forward', i.e. to fix the grid of the real recombination times 
and find all the other saddle point solutions, and the corresponding harmonic number 
A'. The recombination condition Eq. (1.38 i (p^ ii — p' + ip") can be re-written as 
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follows: 

(Ap' + iAp")^ 
Ap' 
Ap" 



2{Nuj~Ip), 

p — ylo sinh(<^r) cos(0r), 



yielding 



{Ap'f - (Ap")^ + 2iAp'Ap" = 2{Noj - Ip 



(1.160) 
(1.161) 
(1.162) 



(1.163) 



Since the right-hand side of this equation is real, we obtain that Ap'Ap" — 0. For 
above threshold hai'monics {Nlj — Ip) > and Ap' / 0, Ap" = 0. For below thresh- 
old harmonics {Nuj — Ip) < and Ap' = 0, Ap" / 0. Separating the imaginary 
and the real parts in the Eqs. ( |1.36|1.37^ , we obtain the four equations quoted in the 
main text, see Eqs. ( |1.43|1.44|1.47|1.48^ . Supplementing these equations for above 
threshold harmonics with Ap" = yields 



P2 ~ smh{(j)r) cos{<l)r) , 
and for below threshold harmonics with Ap' — 0, yielding 

pi — sin((jf>r) cosh{(l)r), 



(1.164) 



(1.165) 



we obtain five equations. 

For above threshold harmonics: for each fixed (j>'r we use the Eqs. ( 1.47|l.48 to 
express cj>i, 4>'i viap2 andpi and then we use Eq. (1.164 1 to exclude p2 and substitute 
</)'^(pi, (/)"), 4>'l{pi,(f)r) and P2(pi, 0") into the Eqs. (1.43 1.44 l Using the gradient 
method, we can now find the minima of the function F = F^ + f| in the plane of 
Pi and <^" for each fixed cf)'^. The minima define the saddle point solutions forpi and 
(f>r. Knowing pi and 0", we find 4)^, <f>'l, p2 from the Eqs. ( |l.47|l.48|l.T64| . Finally, 



the corresponding harmonic number can be calculated from ( Ap') " = 2{Nu! — Ip), 
yielding Nuj = Aq (pi — sin (p'r cosh (f'r)'^ / 2 + Ip. Naturally, the harmonic number 
defined this way does not have to be integer. 

For below threshold harmonics, the procedure is essentially the same. For each 
fixed 4>[., we use the Eqs. ( 1.47 1.48 i to express (f>i, 4>i via p2 and pi and then we 
use Eq. ( 1.165| to exclude pi and substitute 4>iiP2, 0")' 4'i{P2, 4>r) and pi(p2, 4>r) 
into the Eqs. ( |1.43| and |1.44[ l. Using the gradient method, we can now find the 
minima of the function F = Ff -f F| in the plane of p2 and 0" for each fixed 
(f>r. The minima define the saddle point solutions for p2 and (f>r- Knowing p2 and 
(pr we find 0^, (/>'/, pi from the Eqs. ( |l.47|l.48|l.l65| . Finally, the correspond- 
ing harmonic number can be calculated from (Ap")" = 2{Ip — Noj), yielding 
Nuj = Ip - AI (p2 - sinh (/>" cos <^r)^/ 2. 

In this method, it is convenient to determine the return time (f>t7- corresponding to 
the threshold harmonic number Nt = Ip/cj. This can be easily done, since, at the 
threshold, p2 — sinh(<^") cos((;/>5,) and pi — sin(<j!>{,) cosh((;/>"). Thus, we can use 
these equations together with the Eqs. ( 1.47 1.48 ) to express (f>'^, 4>'l via (jij. and 0" 



and then use the Eqs. ( 1.43 1.44 1 to find a minimum of the function F — F( + F2 
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in the plane of ^J, and 0", representing the threshold values of (p'r and <^". Once the 
threshold value 0^^ of <^r is known, one can separately implement the procedures 
described above for below, (b'j. < A'f, and above, S'^ > fif>L, threshold harmonics. 
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